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Abstract 


In geophysical and plasma contexts, zonal flows are well known to arise out of turbu¬ 
lence. We elucidate the transition from statistically homogeneous turbulence without 
zonal flows to statistically inhomogeneous turbulence with steady zonal flows. Start¬ 
ing from the Hasegawa-Mima equation, we employ both the quasilinear approxima¬ 
tion and a statistical average, which retains a great deal of the qualitative behavior 
of the full system. Within the resulting framework known as CE2, we extend recent 
understanding of the symmetry-breaking ‘zonostrophic instability’. Zonostrophic in¬ 
stability can be understood in a very general way as the instability of some turbulent 
background spectrum to a zonally symmetric coherent mode. As a special case, 
the background spectrum can consist of only a single mode. We find that in this 
case the dispersion relation of zonostrophic instability from the CE2 formalism re¬ 
duces exactly to that of the 4-mode truncation of generalized modulational instability. 
We then show that zonal flows constitute pattern formation amid a turbulent bath. 
Zonostrophic instability is an example of a Type I s instability of pattern-forming sys¬ 
tems. The broken symmetry is statistical homogeneity. Near the bifurcation point, 
the slow dynamics of CE2 are governed by a well-known amplitude equation, the 
real Ginzburg-Landau equation. The important features of this amplitude equation, 
and therefore of the CE2 system, are multiple. First, the zonal flow wavelength is 
not unique. In an idealized, infinite system, there is a continuous band of zonal 
flow wavelengths that allow a nonlinear equilibrium. Second, of these wavelengths, 
only those within a smaller subband are stable. Unstable wavelengths must evolve 
to reach a stable wavelength; this process manifests as merging jets. These behav¬ 
iors are shown numerically to hold in the CE2 system, and we calculate a stability 
diagram. The stability diagram is in agreement with direct numerical simulations of 
the quasilinear system. The use of statistically-averaged equations and the pattern 
formation methodology provide a path forward for further systematic investigations 
of zonal flows and their interactions with turbulence. 
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Chapter 1 
Introduction 


Zonal flows are turbulence-driven sheared flows. They are usually associated with a 
direction of symmetry. In planetary atmospheres, they flow along lines of latitude, 
parallel to the equator, and the direction of flow alternates with latitude. In that 
context, zonal flows are associated with the azimuthal symmetry. In magnetically 
confined toroidal plasmas, zonal flows consist of E x B flows produced by toroidally 
and poloidally symmetric fluctuations of electric potential. The direction of flow is 
along a flux surface and varies radially. 

Zonal flows have taken on special significance in plasma physics because they 
are thought to regulate drift-wave turbulence. In particular, evidence is mounting 
that turbulence driven by the ion-temperature-gradient (ITG) instability in toroidal 
plasmas is suppressed by zonal flow or mean shear flow. 1 Furthermore, zonal flows 
are thought to play a role in triggering the L-H transition. The enhanced plasma 
performance of the H-mode is viewed as essential to any viable fusion reactor, and 
zonal flows may play an important part of the H-mode. 

One mechanism by which shear flow is believed to suppress turbulence is shear- 
enhanced decorrelation (Biglari et al. 1990, Diamond et al. 2005, Terry 2000). The 
basic idea is that the flow causes a turbulent eddy to stretch and elongate, making it 
more likely for that eddy to break apart. This reduces the length scale of turbulence 
and hence reduces the resultant turbulent transport. Numerical simulations seem 
to corroborate the idea, directly implicating zonal flows in reducing the levels of 
turbulent fluctuations ( bn et al. 199 ). This simple, powerful idea has been incredibly 
influential, spawning an entire genre of inquiry, and zonal flows have been under 
intense study by the plasma physics community ever since. Any means that might 
help in taming the beast of tokamak turbulence is pursued with vigor. 

Zonal flow is also prominent in geophysical contexts. For example, Figure 1.1 
shows Jupiter with visible zonal bands. The alternating bands flow in alternating 

1 Zonal flow refers to turbulence-driven flow, and it typically oscillates in space with finite radial 
wavenumber. Mean shear flow is caused by diamagnetic effects associated with the mean pressure 
profile. 


1 




Figure 1.1: Jupiter, with zonal bands visible. Image from NASA spacecraft Cassini. 

directions. 2 All of the gas giants in our solar system have zonal flows, not just 
Jupiter. Due to their visibility, the atmospheric science community has studied zonal 
flow for decades. 

Zonal flows and zonal magnetic fields are also beginning to be observed in as- 
trophysical simulations of accretion disc turbulence driven by the magnetorotational 
instability (Johansen et al. 2009, Kunz and Lesur 2013). We cannot currently observe 
and may never be able to directly observe zonal structure of accretion discs, but our 
understanding of their dynamics may hinge upon the behavior of zonal fields. 

Since zonal flows are driven by turbulence, any understanding of zonal flows must 
begin with an understanding of turbulence. In this chapter we start by introducing 
some important aspects of turbulence. Then, we turn to zonal flows and review the 
experimental, numerical, and theoretical literature, separated into geophysical and 
plasma physics sections. 


1.1 Turbulence in Fluids and Plasmas 

The word turbulence conjures up images of chaotic motion, of disorder. Turbulence 
would seem to destroy any semblance of regularity or organization. Typically, smooth, 

2 Animated images are available at http: //ciclops. org/view/92/Jirpiter_Mosaics_and_ 
Movies_-_Rings_Satellites_Atmosphere and http://www.nasa.gov/centers/goddard/ 
multimedia/largest/EduVideoGallery.html. 
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laminar flow such as regular pipe flow or Rayleigh-Benard convection rolls gives way 
to disorder, turbulence, and a jumble of scales. Out of this turbulence, seemingly by 
magic, coherent structures such as zonal flows can form, as we shall see. 

Turbulence theory in fluids and plasmas has varying objectives. In 3D, homoge¬ 
neous, isotropic, incompressible Navier-Stokes turbulence, theory has been trying to 
understand intermittency of the inertial range. In geophysical fluid dynamics, the 
goal of theory is to understand the atmospheres of not only other planets, but also 
our own. The Earth’s combined atmosphere-ocean system constitutes an incredibly 
complex dynamical system, one which determines our climate. In fusion theory, the 
ultimate objective of turbulence theory is to predict and control the level of turbu¬ 
lent transport. When one attempts to build a fusion reactor, out of the many, many 
factors that must be considered, the effect of microturbulence often boils down to a 
single number: the energy confinement time. The greater the level of turbulence, the 
worse the confinement of heat and energy within the plasma. 

In the following sections we introduce a tiny bit of basic turbulence theory. For a 
comprehensive introduction, see Davidson (2004). 


1.1.1 Cascade in 3D 


The natural place to start is with Kolmogorov’s explanation of the energy cascade in 
3D neutral-fluid turbulence ( hisch 1995). The famous Kolmogorov scaling is one of 
the most fundamental and celebrated results of neutral fluid turbulence theory. It is 
one of the first quantitative, successful predictions of fully developed turbulence. The 
result concerns turbulence of the incompressible Navier-Stokes equation, 


— + v ■ Vv = — Vp + ^V 2 v, 
at 

V-v = 0, 


(1.1) 

( 1 . 2 ) 


where p is the pressure divided by density and v is the viscosity. The Kolmogorov 
theory makes a definite prediction for the energy spectrum in wavenumber space. 
First, several assumptions are made: 

1. The turbulence is statistically homogeneous and isotropic. 

2. Energy flows locally in k- space. 

3. There is an inertial range of fc-space where the turbulence does not “know” 
about forcings at the large scale or viscosity at the small scale. 

The average energy density 3 E and omnidirectional energy spectrum E(k) of the 
flow are related by 



(1.3) 


3 Using the density instead of the total energy prevents us from having to deal with infinite 
energies in an infinite fluid. 
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The energy density E is decomposed as a sum over wavenumbers. The spectrum 
E(k) depends only on the magnitude of the wavenumber, k = |k| (where = denotes 
a definition), due to the isotropy assumption. 

Physically, one often speaks of “eddies” in a turbulent flow. The typical picture 
is that energy is somehow injected into the system at large scales, perhaps due to 
mechanical stirring of the fluid, and gives rise to eddies. These turbulent eddies 
interact with each other in some way, giving rise to smaller scale eddies. Energy 
flows from the larger scales to the smaller scales in this scenario. This is related to 
assumption 2. Eventually, when energy reaches small enough scales, viscosity becomes 
important and the energy is dissipated. The physical attributes at an intermediate 
scale are assumed to not depend on the precise behavior at the large or the small 
scales. 

Let us make this more precise. Energy is injected at a rate £ at the large scales, 
the forcing scales, designated by wavenumber kf. Energy is assumed to flow locally 
through fc-space without dissipation until it reaches the viscosity-dominated small 
scales, designated by wavenumber k v . In a statistically steady state, the energy flux 
through every scale must, on the average, be £, until it is dissipated. At intermediate 
wavenumbers, kf <C k <C k u , the locality assumption means that the turbulence 
cannot depend upon kf or k ui but only on the scale k and the energy flux e. There 
are no other local quantities it can depend on. 

These ideas may be expressed as an advection equation in k- space. 4 While the 
Kolmogorov argument is essentially dimensional and not a quantitative calculation, 
the advection equation is handy for systematizing the assumptions and tracking the 
dimensions of all quantities. The assumptions lead one to being able to write the 
advection equation 


dE(k,t ) d (Ak , A ... . . . 

+ * ( aA M) ) = rf( 7) ' (1 ' 4) 

And in a statistically steady state, E(k,t ) will not depend on t. This is a local 
conservation equation, where forcing but not dissipation has been built in. This 
equation can be considered as part of a Fokker-Planck equation, where Ak/At is the 
“drift velocity” through k- space. We will consider this equation at some k larger than 
kf. 

When speaking of a given scale k or l ~ A; -1 , it will be convenient to assign a 
width to the scale. It is most natural to break the scales up logarithmically. That is, 
starting from the largest scale kf 

k f —> 2 1 k f —» 2 2 k f —» 2 3 k f —> 2% —> • • • . (1.5) 

A given scale labeled A’ can be considered to contain wavenumbers from k/2 to k. 
By the locality assumption, the Ak appearing in the advection equation can only be 
k. 

4 The author first learned of this approach from G. Hammett. 
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At is what is called a nonlinear correlation time, or an “eddy turnover time” r k . 
An eddy turnover time is defined as the time it takes for a fluid element of speed v k 
to cross an eddy of size l ~ /c -1 , 

r k = —. ( 1 . 6 ) 

Vk 

Here v k is the characteristic speed of eddies of size / and is quantified through 


>1= f dk'E(k') 
Jk/2 

~ kE(k), 


(1.7) 

( 1 . 8 ) 


which gives v k ~ \JkE{k) (ignoring constants of order unity). 

Now, integrate the advection equation from k = 0 to some k > kf. The energy 
spectrum E(k) is assumed to vanish at k — 0. One finds 


—E(k)=s, (1.9) 

k 2 viE(k ) = e, (1-10) 

k 5/2 E 3/2 = £ . (i.n) 

One thus obtains the Kolmogorov scaling for the inertial range, 

E(k) = C£ 2/3 k~ 5/3 , (1.12) 

where C is simply an order-unity constant. 


Kolmogorov Scaling from Pure Dimensional Analysis 

Another way to obtain the Kolmogorov scaling is through dimensional analysis with¬ 
out any recourse to the physics. This approach yields less intuition than the physical 
picture of eddy turnover, but it is a useful demonstration of the power of dimensional 
analysis. The locality hypothesis demands 

E(k) = f(£,k). (1.13) 

The dimensions of the average energy density, energy spectrum, and energy flux 
are given below, where L is the dimension of length and T is the dimension of time. 

(1.14) 

(1.15) 

(1.16) 


~ L 2 

E ~ X2> 

E ( k ) ~ 

L 2 

£ rs - / - 

T 3 
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log T 



Figure 1.2: Time scale of inertial effects (blue, solid) and viscous effects (red, dashed) 
as a function of wavenumber. 


Now, suppose the function / ~ £ a k°. Then the dimension of E would be 

„ L 3 ( L 2 \“ 1 

K [ — ) —. 

T 2 V T3 ) L/3 


(1.17) 


Satisfying dimensional consistency requires a = 2/3 and f3 = —5/3. Thus, E ~ 
e 2 / 3 A; -5 / 3 is recovered. 


Energy Dissipation and the Viscous Scale 

At a given k in the inertial range, the eddy turnover time r k is given by 


Tk ~ 


1 

kv k 


k~ z/2 E~ l/2 ~ e -1 / 3 k~ 2/3 . 


(1.18) 


From the Navier-Stokes equation, the timescale for viscous processes at a scale k can 
be seen to be 



(1.19) 


The process with the shorter time scale dominates. Dissipation becomes important 
at the scale Ay where = t\.. For k < Ay, inertial effects dominate, while for k > Ay, 
viscous effects dominate (see Figure 1.2). 

Setting = r k gives an estimate for the viscous scale, or Kolmogorov scale: 


/ £ \ 1/4 /V \ 1/4 

Ay=(^-J or . (1.20) 

The dissipation rate is E — (w ■ V 2 v). We have already assumed that it acts 
only at scales smaller than or comparable to the Kolmogorov scale l v . If we look at 
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the Kolmogorov scale, then substituting V ~ k u , we hnd 



This shows that the dissipation acts primarily at the Kolmogorov scale; nothing much 
is happening at smaller scales. This result is also important because it shows that 
dissipation is independent of the viscosity, even as v —» 0. 


1.1.2 Cascade in 2D 


The nature of cascades are different in two dimensions (Kraichnan 196' ). Instead 
of just the energy being conserved, in 2D there are two quadratic quantities that 
are conserved by the nonlinear interaction: energy and enstrophy. The 2D case is 
important because geophysical flows are quasi-2D due to atmospheric stratification 
(Pedlosky 1987, Vallis 2006), and plasma flows are quasi-2D due to the magnetic held. 

Again assume statistical isotropy and homogeneity. Instead of forcing at large 
scales as in 3D, assume that forcing occurs at some intermediate length scale or 
wavenumber. Then there is a dual cascade, with two inertial ranges rather than 
just one. Energy cascades from the forcing scale to larger scales, whereas enstrophy 
cascades from the forcing scale to smaller scales. The energy cascade is called the 
inverse cascade, while the enstrophy cascade is called the direct cascade. The energy 
spectrum in the inverse cascade range is E{k) = £ 2 / 3 fc~ 5//3 , where e is the energy 
flux through wavenumber space. In the direct cascade range, the energy spectrum is 
E(k) = rj 2 ^ 3 k~ 3 , where r/ is the enstrophy flux through wavenumber space. 5 

The flow of energy to large scales and the flow of enstrophy to small scales can 
be understood as a consequence of the conservation laws (Fjprtoft 1953, Kraichnan 
1967). There is some energy spectrum, E{k), with total energy density given by 



( 1 . 21 ) 


Let Z{k) be the enstrophy spectrum. It is related to the energy spectrum by Z{k) = 
k 2 E{k), so that the total enstrophy density is 



( 1 . 22 ) 


In other words, the enstrophy is weighted by a higher power of wavenumber than 
the energy is. Vallis ( 006) showed that if the energy spectrum spreads out under 
the constraint of conservation of both total energy and enstrophy, then the centroid 


of the energy spectrum must move to smaller k (larger scales) and the centroid of 


5 Kraichnan (1971) showed that a logarithmic correction needs to be applied in the enstrophy 
inertial range; see also Bowman (1996). 
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the enstrophy spectrum must move to larger k (smaller scales). The tendency for 
energy to accumulate at large scales will be especially important for understanding 
the generation of zonal flows in geophysical contexts. 

1.1.3 Statistical Theories of Turbulence 

The statistical approach to understanding turbulence, which this thesis takes, com¬ 
plements other methods such as making detailed measurements of plasma fluctuations 
or performing direct numerical simulations (DNS). Those methods can accumulate 
reams of data so vast that it can be unclear how one should go about making sense of 
it all. The aim of the statistical approach is to focus on the macroscopic quantities of 
interest, such as transport coefficients, energy spectra, and the like. By working with 
averaged quantities from the outset, one can circumvent the rapid spatiotemporal 
fluctuations and potentially see a clearer view of the physics. Of course, there is no 
free lunch. As a consequence of averaging a nonlinear equation, one is generally left 
with the average of an unknown quantity: a closure problem. Various statistical clo¬ 
sures, perhaps motivated by physical considerations, provide different approximations 
for the unknown terms. A major difficulty with this approach is that the closures 
are essentially uncontrolled approximations; the nonlinearity inherent to turbulence 
makes it hard to know exactly what is lost. The closure might obliterate some highly 
coherent or correlated phenomena. Nevertheless, these difficulties do not invalidate 
the statistical approach, from which much has been learned (Frisch 1995, Ivraich- 
nan 1959, 1964b, Krommes 2002). Historically, the majority of theoretical studies 
into turbulence that follow this approach assume homogeneous statistics, where the 
statistics of turbulent quantities do not depend on position. Consequently, most of 
the theoretical machinery that has been developed also applies only to homogeneous 
statistics, with comparatively little effort devoted to inhomogeneous statistics. The 
main line of work in this thesis involves inhomogeneous statistics. 

1.2 Zonal Flows 

The simplest model in which zonal flows arise naturally out 
system 

d t w + v • X7w + — f + D, 

where / is a forcing term, D represents dissipation, and 

w = VV (1-24) 

Here, ifj is the stream function, v = z x is the velocity, and w = z ■ V x v is the 
vorticity. This equation will be discussed much more fully in Chapter 2. The equation 
is often used as the simplest, most reduced description of atmospheric turbulence. The 
behavior of turbulence and zonal flow even in this simple system is still studied today. 
In this introductory chapter, we use this equation to highlight a few key points. 


of turbulence is the 2D 
(1.23) 


1.2.1 Zonal Flows in Geophysics 

We briefly 6 review some of what is known about zonal flows in geophysical contexts. 
For more information, see the works of Pedlosky (1987), Vallis (2006), Vasavada and 
Showman (200 ) and references therein. 

Jupiter, for example, has prominent, easily visible zonal jets. It has roughly 30 
zonal jets, and they have been remarkably stable over time. Measurements by Voyager 
in 1979 and Cassini in 2000 indicate the zonal wind profile has barely changed in that 
time period. Compared to Jupiter’s equatorial radius of 70,000 km, we can directly 
observe at most only a few hundred kilometers into the atmosphere. Little is known 
about the turbulence and zonal wind deeper down. In the upper atmosphere, zonal 
jet speed is mainly measured by assuming that clouds are passive tracers of the 
zonal wind. This is not perfect, due to for example, larger clouds averaging over 
an extended spatial region, but it seems to be somewhat successful. This technique 
can only measure jet speed at cloud level. The energy source of the zonal jets is 
hypothesized to be buoyant convection from a hot planetary interior ( /asavada and 
Showman 2001 ). On Earth, zonal flows occur can occur in both the ocean and 
atmosphere, but the flows tend to meander with complex dynamics, and there are 
not as many jets. 

One idea deserves special note. The notion of the Rhines scale has been enor¬ 
mously influential in the geophysical literature ( Urines 1975). The Rhines scale pur¬ 
portedly estimates the jet width or spacing, and is given by 



where U is the rnrs velocity and (3 is the northward gradient of the Coriolis parameter. 
Inversely, we can express the characteristic Rhines wavenumber as 

k R = {P/U) 1 ' 2 . (1.26) 

We give a couple of ways of obtaining the Rhines scale (Vallis and Maltrud 1993, 
Vasavada and Showman 2005). The first method is essentially dimensional analysis. 
Let 

v = nx + vy (1-27) 

and w = z ■ V x v = d x v — d y u. Then (1.23), rewritten here as 

d t w + v ■ Vtc + v/3 = f + D, (1-28) 

can be used to find the Rhines scale by heuristically balancing the magnitudes of 
the Rossby wave term (the /3 term) and the nonlinear advection term. If we treat 
u ~ v ~ U and k x ~ k y ~ k, then we find that the nonlinear advection term is 
roughly k' 2 U 2 , and the linear term is roughly Uf3. Where these are equal gives this 
Rhines scale. 

6 Very briefly, since this is not the author’s area of expertise. 
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A slightly more refined analysis would allow for the zonal jets to have a different 
magnitude and length scale than the eddies. Let v and ( be the zonally-averaged 
velocity and vorticity, U be the characteristic velocity of zonal flow, and u',v 7 be the 
characteristic velocity of the eddies. We suppose U u',v'. Then 

v • VC ~ v 7 ■ VC + v ■ VC' (1.29) 

= -v'dyu(y) + u(y)d x (d x v' - d y u ,'). (1.30) 

If we look particularly at the first term of (1.30),' then we see that the advection term 
goes like vk^U, which can be compared with v/3. Equating them gives the Rhines 
scale. 

A more physical argument views the Rhines scale as a transition scale between 
the regimes where inertial, isotropic turbulence and Rossby-wave activity dominates 
(Rhines 1975, Vasavada and Showman 2005). Assume the turbulence is forced at 
small scales. At wavenumbers greater than the Rhines scale, the eddy-turnover time 
scale is shorter than the time scale of Rossby waves, so standard 2D turbulence results 
with an inverse cascade. Energy proceeds towards larger scales until it reaches the 
Rhines scale. Then the time scale of Rossby waves becomes shorter than the eddy 
turnover time, so Rossby waves dominate. The idea is that these large scale waves 
are inefficiently forced by the turbulence, and energy cannot easily cascade to length 
scales larger than the Rhines scale, so energy piles up at k r and the inverse cascade 
slows. The turbulent frequency is roughly Uk and the Rossby wave frequency is 
oj = —k x /3/k 2 . Let (p be the angle between east (the x direction) and the direction of 
wave phase propagation, so that cos 0 = k x /k. Then equating the turbulent frequency 
and the wave frequency leads to an anisotropic Rhines scale, 

k 2 R — Jj \cos 0|- (1-31) 

A plot of this anisotropic “dumbbell” shape is shown in Figure 1.3. The dumbbell 
outline is where the inverse cascade halts. The anisotropic shape offers an explanation 
for why energy piles up on the k y axis where k x = 0, leading to the preference of 
zonally-symmetric structures. This scenario appears to have been confirmed (Vallis 
and Maltrud 1993), although some have called it into question by arguing that the 
small scales directly force the zonal flows (Huang and Robinson 1998). 

Large-scale friction or drag is critical for getting the Jovian jets correct ( 7 asavada 
and Showman 2005). Without it, energy would slowly leak past the Rhines scale 
into larger scales. If this energy is not damped somehow, then on a long enough 
time scale, even the large scales would isotropize and distinct zonal jets would not 
exist. Large-scale friction damps the energy that would leak past the Rhines scale. 
However, if friction were too large, energy would damp before it could cascade up to 
the Rhines scale, and so no jets would form and isotropic turbulence would result. 

7 This is not particularly justified without further argumentation. Based on magnitudes, one 
might expect the other term to dominate because the length scale of turbulence is smaller and so 
its spatial derivatives are larger. 
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Figure 1.3: Anisotropic Rhines scale. Outside the dumbbell, at high k, inertial tur¬ 
bulence has a shorter time scale and dominates. Inside the dumbbell, Rossby wave 
dynamics are faster. The hypothesis is that the inverse cascade cannot penetrate into 
the dumbbell, so energy piles up on the k y axis where k x = 0. 


1.2.2 Zonal Flows in Plasmas 

Our definition of zonal flow in plasma refers only to the zero-frequency flows. We 
exclude geodesic acoustic modes (GAMs) ( Vinsor et al. 1968) from our definition of 
zonal flow. Some authors describe GAMs as oscillatory zonal flows, but in this thesis 
we do not. 


Theory and Simulations 

Efforts at developing a systematic theory of zonal flows have almost exclusively fo¬ 
cused on the scenario where the zonal flows are assumed to be long wavelength com¬ 
pared to the scale of the turbulence (Connaughton et al. 2011, Diamond et al. 2005, 
Krommes and Kim 2000, Smolyakov et al. 2000b). This remains true despite the 
fact that in simulations and experiments, zonal flows tend to be of scale compa¬ 
rable to that of the turbulence. When the long-wavelength zonal flow assumption 
is made, the resulting interaction of turbulence and zonal flow can be described in 
terms of a wave kinetic equation. In this type of description, one might imagine 
a sea of drift-wave packets evolving in an weakly-inhomogeneous medium of zonal 
flows. The wave-kinetic formulation has intuitive advantages because the turbulent 
wave action is materially conserved along phase-space trajectories. In addition to 


11 







the long-wavelength assumption, many studies make a single-harmonic assumption 
where only a single Fourier mode of the zonal flows is retained (Connaughton et al. 
201 1). When that is done the zonal flow wavelength cannot be found from the theory 
but is left as an undetermined parameter. The state of things indicates that theory 
of zonal flows is still in its infancy. 

Some studies have focused on a generalized modulational instability, in both the 
geophysics and plasma literature (Connaughton et al. 2010, Gill 1974, Lorenz 1972, 
Manin and Nazarenko 1994, Smolyakov et al. 2000b, Wordsworth 2009). 8 In analytic 
studies of these instabilities, typically a single eigenmode, referred to as the primary 
wave, is used as the background upon which the perturbation grows. For example, in 
a periodic box, a single Fourier mode is an exact solution to the nonlinear vorticity 
equation. A conceptually close cousin of modulational instability is secondary insta¬ 
bility (Plunk 2007, Pueschel et al. 2013, Rogers et al. 2000). In secondary instability, 
a growing linear eigenmode, the primary mode, acts as a background upon which a 
secondary perturbation grows. If the secondary mode grows much faster than the 
primary, then the primary can be treated as stationary. 

Other simulations have investigated various aspects of zonal flows. For instance, 
Nakata et al. (2012) examined entropy transfer via zonal flows in gyrokinetic sim¬ 
ulations. Along the same lines, Makwana et al. (201 1, 2012) looked at how zonal 
flows interact with damped modes to regulate turbulence. Xanthopoulos et al. (2011) 
studied the effects of the magnetic equilibrium in stellarator geometry. Waltz and 
Hollanc (2008), after turning off drift wave-drift wave nonlinear couplings, concluded 
that the drift wave-zonal flow coupling accounts for most of the nonlinear satura¬ 
tion. Other models, using fluid equations and simplified geometry, are also a fruitful 
ground with which to gain intuition and insight. The Hasegawa-Wakatani system is 
one such model which has been used to study zonal flows (Hasegawa and Wakatani 
1983, 1987, Pushkarev et al. 2013). Numata et al. (2007) first studied the Modified 
Hasegawa-Wakatani system, which corrects the treatment of zonal flows when the 
equations are restricted to two dimensions. 

Experimental Observations 

It is tough to make direct measurements of zonal flows in plasmas. First, the zonal flow 
involves only an electric potential and plasma flow, which are not easily observed. In 
contrast, the GAM is associated with an rn — 1 poloidal density fluctuation which can 
be measured more readily. Furthermore, the GAM oscillates at moderate frequency, 
whereas the zonal flow fluctuates at zero or low frequency, which is also difficult to 
measure by certain techniques. Zonal flows are often not zero frequency in practice, 
but fluctuate on a much slower time scale (a few kHz) than the turbulence (tens of 
kHz). 

The main diagnostic tools used to measure zonal flows are the Langmuir probe, 
the heavy ion beam probe (HIBP), beam emission spectroscopy (BES), and Doppler 
reflectometry (Estrada, Fujisawa 2009). 

8 Generalized in the sense that it is not restricted to the original meaning of long-wavelength 
modulations that vary in the same direction as the primary wave. 
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Langmuir probes measure the ion saturation current and floating potential 
(Hutchinson 200"). Fluctuations in floating potential are usually analyzed as fluctua¬ 
tions in plasma potential and fluctuations in ion saturation current as fluctuations in 
plasma density. Potential measurements in multiple locations can be used to calculate 
electric fields and hence ExB flows. One signature of zonal flow is correlation in 
electric potential between positions on the same flux surface but separated toroidally 
or poloidally. Density measurements are useful because finding a weak correlation 
in density fluctuations while detecting a strong correlation in potential fluctuations 
enhances one’s confidence that the observed phenomenon is in fact a zonal flow. 
Langmuir probes are restricted to cooler plasmas because the probes would otherwise 
not survive, so Langmuir probes cannot be used for measurements in the core of 
high-performance plasmas. Owing to their simplicity, Langmuir probes are widely 
used when feasible. 

The HIBP diagnostic provides a direct measurement of plasma potential, even 
in the plasma core (Crowley 1994, Ido et al. 2002). Heavy, singly charged ions are 
injected into the plasma at high energies (hundreds of keV). Upon impact with elec¬ 
trons, some ions undergo ionization into a double charge state, and these so-called 
secondary ions are deflected more strongly in the magnetic field. The secondary ions 
also gain energy at the ionization point due to the increase in potential energy from a 
higher charge state. The plasma potential 0 at the ionization point can be determined 
by measuring the difference in kinetic energy between primary and secondary ions at 
a detector. 

The BES diagnostic measures local density fluctuations ( ? onck et al. 1990). A 
neutral beam is injected into the plasma and undergoes collisional fluorescence. The 
emitted light is approximately proportional to the local density. The Doppler shift of 
the emitted light due to the beam’s velocity allows for the separation of the beam H a 
or D a emission from the bulk plasma emission. High spatial and temporal resolution 
is possible with a 2D imaging system. Flow velocity can be calculated from the 
motion of turbulence structures between poloidally separated channels using time- 
delay-estimation techniques. 

Doppler reflectometry yields measurements of both flow and density fluctuations 
( lirsch et al. 2001). Unlike traditional reflectometry, Doppler reflectometry uses an 
angle between the incoming microwave beam and the cutoff layer. Scanning the 
tilt angle allows wavenumber-resolved turbulence measurements. Flow velocity is 
measured from the Doppler shift of the scattered signal. This diagnostic can provide 
high temporal and spatial resolution. 

The first direct observation of zonal flows in the core region of a toroidal plasma 
came from a dual-HIBP measurement in the CHS stellarator ( ujisawa et al. 200 z ). 
Zonal flows were then found in the core of a tokamak plasma in DIII-D with BES 
(Gupta et al. 2006). In both of these cases, the measured zonal flows exhibited a short 
radial length scale comparable to that of the turbulence. Many other identifications of 
zonal flows can be found in the reviews of Fujisawa (2009) and Cstrada and references 
therein. 
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L H Transition 

The transition from L-mode to H-mode has been the subject of intense interest since 
its discovery ( .Vagner et al. 1982). H-mode is associated with a transport barrier 
and a reduced level of turbulence, along with a sharper plasma pressure gradient 
which improves performance. A number of studies have implicated sheared E x B 
flows in the L-H transition, although there is not yet a detailed understanding of the 
associated physics. 

Measurements of E r find that the E x B flow varies rapidly during the L-H 
transition, whereas the pressure profiles and the resulting diamagnetic flow take longer 
to evolve. It also found that the increase in E r shear occurs before the decrease 
in turbulent fluctuations, consistent with shear flow causing turbulent suppression 
(Burrell 1999, Estrada et al. 2009, Meyer et al. 2011, Moyer et al. 1995). 

With recently improved spatiotemporal resolution, many devices have observed 
between L-mode and H-mode an intermediate, transient phase, which is called I-phase 
(Colchin et al. 2002, Estrada et al. 2012, 2010, 2011, Schmitz et al. 2012, Xu et al. 
2011). 9 The I-phase is characterized by oscillations in the zonal flow and turbulent 
fluctuations. These oscillations often show a characteristic predator-prey behavior, 
with the E x B flow (the predator) following the density fluctuations (the prey) with 
a phase delay of 90°. 

On the theoretical side, no first-principles simulation has reproduced the L-H 
transition. Consequently, theoretical investigations have focused primarily on reduced 
models with various assumptions and approximations. Initially, these studies were 0D 
and modeled the predator-prey interaction only between the mean shear flow and the 
fluctuation level (Diamond et al. 1994). Later, Kim and Diamond (2003) extended 
the model to include suppression of turbulence by both mean flow and zonal flow. 
This two-predator, one-prey model exhibits pre-transition oscillations. In the model, 
the zonal flow triggers the transition, and then the steep-gradient-driven mean flow 
sustains the H-mode. That work has been developed further into a ID, radially- 
extended model that evolves turbulence intensity, zonal and mean flow shear, and 
pressure and density profiles ( : iki et al. 201 ). 

1.3 Overview of this Thesis 

The strategy of this thesis is to start at the basics and develop a systematic theory 
of zonal flows from the bottom up. To this end, we use the simplest models in order 
to develop a sound theoretical foundation. In Chapter 2, we introduce the Charney- 
Hasegawa-Mima equation, which serves as the model for almost all of the work in 
this thesis. This equation neglects many of the realistic effects in plasmas and fluids, 
which allows for tractable analysis. We also describe the quasilincar approximation 
and the CE2 statistical framework. The quasilinear approximation denotes that the 
fields of interest are divided into a mean field and an eddy field, and then in the 

9 The I-phase should not be confused with the I-mode regime first observed on the Alcator C-Mod 
tokamak (Whyte et al. 2010). 
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equation for the eddy field the eddy-eddy nonlinearity is neglected. We perform all of 
our analysis within the context of the quasilinear approximation. While it is clearly 
not a realistic approximation in all cases, numerical simulation provides convincing 
support that the qualitative behavior, at least of zonal flows, is similar as in the full 
model. This approximation leads naturally to the CE2 statistical framework. Because 
statistical theories of turbulence average over small-scale fluctuations, these theories 
have the advantageous feature of allowing for a steady-state, statistical description 
of a turbulent equilibrium. For example, many theories, including CE2, describe 
turbulence in terms of a two-point correlation function. 

Chapter 3 contains the main physics content of the thesis. We first review the 
recently discovered zonostrophic instability (ZI). In ZI, a statistically homogeneous 
turbulent state that is on average uniform in space becomes unstable to a inhomoge¬ 
neous perturbation. These perturbations grow into saturated zonal flows. We draw 
connections between ZI and modulational instability. Then we show that zonal flow 
can be interpreted as pattern formation, and we expand upon the insights that brings. 
For instance, as a control parameter is varied and the homogeneous turbulent state 
becomes ZI unstable and a new stable inhomogeneous state appears, the bifurcation is 
described by a simple equation with universal behavior. One immediate consequence, 
previously remarked in scattered observations but never explicitly understood math¬ 
ematically, is the existence of multiple solutions to the CE2 equations with varying 
ZF wavelengths. In other words, the width of the zonal jets is not unique; we derive 
this mathematically. We analytically calculate the bifurcation at which zonal flows 
appear and verify it numerically. In Chapter 4, we solve CE2 numerically to find 
equilibria of nonlinearly interacting turbulence and zonal flows. To do this, we use 
Newton’s method to directly solve the steady-state CE2 equations. This technique is 
common for pattern-forming systems. Once the equilibria are found, we also calculate 
their stability. In terms of the ZF wavelength, the region of stability calculated from 
CE2 is consistent with the results of direct numerical simulation of the QL equations. 

Chapter 5 proposes a simple closure for homogeneous turbulence. We are able to 
examine in detail the stability of solutions to the closure, a property mostly neglected 
in the literature. As it stands, this chapter is somewhat separate from the rest of the 
thesis. But if extended, it could be connected back to zonal flows and provide a way 
for further analytic progress. 

Finally, in Chapter 6, we explore directions for future research. For instance, 
we discuss more realistic turbulence closures than the quasilinear approximation and 
CE2, such as the DIA. Such closures are needed to account in some way for the eddy- 
eddy nonlinearities. These more sophisticated approaches would allow for better 
quantitative and qualitative accuracy. Additionally, in this thesis we dealt with the 
Charney-Hasegawa-Mima equation where turbulence is forced by external drive, but 
most models of plasma turbulence relevant to fusion have an intrinsic instability. One 
simple way to proceed along this path is to extend the closure described in Chapter 5 
to allow for inhomogeneity. Then one could perform a similar bifurcation analysis to 
that in Chapter 3. The results of this thesis provide a theoretical foundation for those 
more sophisticated models, but research is needed to understand those situations in 
detail. Toroidal geometry presents a challenge, especially for analytic work. One 
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Tabic 1.1: Important symbols, their meaning, and the equation in which they are 
first used. 


Symbol 

Meaning 

Equation 

= 

Definition 


V 

Viscosity 

(1.1) 

w 

Generalized vorticity 

(1.23) 


Stream function 

(1.23) 

P 

Planetary vorticity gradient (or plasma density gradient) 

(1.23) 

L d 

Deformation radius (or plasma sound radius) 

(2.3) 

OLZF 

Modifies relation of w and (is either 1 or 0) 

(2.6) 

h 

Friction 

(2.8) 

h 

Hypervisocity factor 

(2.8) 

7 

Fundamental dimensionless paramater 

(2.13) 

U 

Zonal flow velocity 

(2.17) 

W 

Covariance of vorticity 

(2.20) 

T 

Covariance of stream function 

(2.21) 

F 

Covariance of random, external forcing 

(2.21) 

x,y 

Difference coordinates of 2-point correlation function 

(2.21) 

y 

Sum coordinate of 2-point correlation function 

(2.21) 

u± 

U{y±\y) 

(2.21) 

V 2 

V 2 - Lf = dl + - Lf 

(2.21) 

/ 

1 — azFL d 2 d ¥ 2 

(2.21) 

t 

e + l- 2 

(3.7) 

Q 

Wavenumber of zonal flow 

(3.8) 

X 

Eigenvalue (i.e., growth rate) 

(3.8) 

t 

q 2 + azpL d 2 

(3.11) 

mnp 

Fourier coefficients of W(x,y y) 

(4.2) 

u p 

Fourier coefficients of U ( y ) 

(4.2) 


plausible path forward is to use CE2 to describe both ZFs and geodesic acoustic 
modes (GAMs) together. Just as we have gained definite insights into the behavior 
of ZFs, we are optimistic that similar insights are possible for GAMs. 

1.4 Mathematical Conventions 

A table of important mathematical symbols is given in Table 1.4. 

1.4.1 Coordinate Convention 

The geophysical communities and plasma communities use opposite coordinate con¬ 
ventions for the direction of inhomogeneity in two-dimensional planar models. In 
planetary atmospheres, zonal flows run in the east-west (x) direction and vary in 
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Figure 1.4: Local coordinate system on a sphere. 


the north-south (y) direction (see Figure 1.4). In tokamaks, zonal flows run in the 
poloidal direction. If one imagines a small box placed at the outboard midplane of 
a tokamak, the poloidal direction becomes the y direction and the radial direction 
becomes the x direction (see Figure 1.5). These opposite conventions for the direction 
of zonal flow require us to make a choice as to which will be followed. We follow the 
convention used in the geophysical community. We do this because the method of 
approach used in this thesis is closely related to recent works in the geophysical litera¬ 
ture. It is significantly easier to comprehend the literature when the same convention 
has been used everywhere. Thinking in terms of the usual tokamak convention then 
requires flipping only a single mental switch. If the alternative of using opposite 
conventions had been chosen, then the practitioner must separately assimilate each 
equation that is encountered. However, as a compromise, we also give the rule to 
transform between conventions for key equations. As an aside, it might also be noted 
that the tokamak convention could be made consistent with the geophysical conven¬ 
tion, if one were to place the region of interest not at the outboard midplane but at 
a poloidal angle of 90°. 

1.4.2 Fourier Transform Convention 

Throughout, we use the convention 


f(k) = 1 dxe-' k *f(x) 

(1.32) 

/M = ^/ dke‘ k *f(k) 

(1.33) 


Often, we work only in Fourier space and drop the hat. 
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Figure 1.5: Local coordinate system at the outboard midplane of a tokamak. (Image 
from the GYRO team.) 
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Chapter 2 

Equations of Motion 


2.1 Modified Hasegawa—Mima Equation 

We take as the starting point the Modified Hasegawa-Mima equation ( homines and 
Kim 2000, Smolyakov et al. 2000a). The Hasegawa-Mima equation (HME) has been 
studied for decades as a paradigm of electrostatic turbulence (Hasegawa and Mima 
1978, Horton and Hasegawa 1994). The basic physics involve a nonuniform back¬ 
ground density profile and the motion of charged fluid elements due to E x B and 
polarization drifts where the fluctuating electric field is self-consistently determined. 
The HME was originally derived from a fluid perspective using the Braginskii equa¬ 
tions (Hasegawa and Mima 1978). It can also be derived in a simple way using a cold 
ion limit of the gyrokinetic equation ( Homines 2006). 

The Modified Hasegawa-Mima (mHME) equation fixes a defect of the original 
version in its treatment of zonal flows. Built into the HME is the assumption of 
adiabatic electrons. But the adiabatic electron response relies on the fast motion 
of electrons along magnetic field lines. Electrons can rapidly flow along magnetic 
field lines to neutralize charge imbalance, but cannot flow across field lines in the 
same way. For fluctuations that are constant on a magnetic surface (i.e., with kn = 
0), the adiabatic electron response breaks down. Therefore zonal flows, which by 
definition are produced by electrostatic fluctuations constant on a magnetic surface, 
are not treated correctly in the HME (Hammett et al. 1993). The mHME modifies 
the electron response to be more physically correct. 

In a 2D formulation, the mHME is typically written as 

d t w(x, y) + v ■ Viu — ndycj) — f + D, (2.1) 

where x corresponds to a radial-like direction, y to a poloidal-likc direction, (j) = 
(L n /p s )e(p/T e is the normalized electrostatic potential, L n is the density gradient 
scale length, p s is the sound radius, T e is the electron temperature, w = V 2 0 — acj) is 
the generalized vorticity and is related to ion gyrocenter density fluctuations 5nf by 
w = — (L n /p s )5nf /no where n 0 is the background density, a is an operator that is zero 
when acting on zonal flows and unity when acting on drift waves, the magnetic field 
is in the z direction, v = z x V</> is the ExB velocity, n is related to the density scale 
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length, / is some kind of forcing that drives turbulence, and D is a dissipation term. 
Lengths are normalized to the sound radius p s and times are normalized to the drift 
wave period u j ~ 1 = (L n /p s )Q,~ l . These normalizations and scalings are convenient to 
make w, 0, and the active length and time scales of order unity. Additionally, they 
also allow us to^set n = 1. 

The terms / and D produce forced, dissipative turbulence. Many studies examine 
the ideal limit in which both / and D are neglected ( lorton and Hasegawa 1994). The 
systems in such cases are Hamiltonian and conserve an infinite number of quantities. 
The difference between the ideal limit and the non-ideal limit here is much the same 
as the difference between the neutral-fluid Euler equation and Navier-Stokes equation. 
Studies of the ideal limit may yield qualitative insight regarding statistical equilibrium 
or cascades of conserved quantities ( ee 1952, Zhu and Hammett 2010) or may reveal 
the tendency of a system to form coherent structures. In contrast, this thesis is 
concerned with forced, dissipative turbulence. 

Since a fixed, spatially-independent profile gradient is used, the model is called 
a local model, as described in Chapter 1. One imagines the domain of the system 
is a small box within the much larger physical system. The local approach tries 
to extract as much physics as possible from as simple a system as possible. The 
local approach can often be justified in terms of the smallness of the ratio p/L , 
where p is the gyroradius of the relevant species, and L is the length scale of the 
macroscopic parameter, e.g., for density, L” 1 = \d\nno/dx\. This approach retains 
the physics of the existence of a gradient in density (or other macroscopic parameter) 
but does not require detailed spatial profile information. This approach can capture 
a great deal of the physics involved and also remove the necessity of dealing with 
complicated boundary conditions. Since turbulence has a small length scale, the 
usual argument goes that regardless of the boundary conditions that are used in 
theory or simulation of the model equations, a few correlation lengths away from the 
boundaries the turbulence should not be affected by the boundaries. For simplicity, 
many simulations use periodic boundary conditions. Even when the local approach 
cannot be rigorously justified in an asymptotic ordering for some realistic situation, 
it is still a useful method for obtaining qualitative insight. 

2.2 (Equivalent) Barotropic Vorticity Equation 

It has been long known that the HME is mathematically very similar to an equation 
that arises in a geophysical context (Pedlosky 198i ). When one writes the equation 
of motion for an incompressible fluid on the surface of a rotating sphere, one has 
what is known as the quasigeostrophic (QG) equation for barotropic vorticity. This 
2D formulation on the surface of the sphere is a decent qualitative approximation, 
because due to the rotation the oceans and atmosphere are stratified into horizontal 
layers. In an approach which is in very much the same spirit as the local approach 
described above, theoretical geophysicists consider not a rotating sphere but a /3 
plane. A /3-plane approximation simplifies the model by linearizing the variation 
of the Coriolis term in the equation of motion. In the full spherical geometry, the 
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Coriolis term varies nonlinearly (that is, sinusoidally) over the latitude of the sphere. 
The /3-plane approximation retains the fact that there is variation, but keeps only a 
linear variation. With such an approach, Rossby waves (the analog of plasma drift 
waves) are easily analyzed. Sometimes the geophysical literature keeps in the /3-plane 
approximation something known as the deformation radius L d , which is the length 
scale at which rotational effects become as important as buoyancy or gravity waves. 
Even though L d is usually thought of as being comparable to turbulent length scales, 
many theoretical studies continue to neglect gravity-wave effects by taking infinite 
L d (Scott and Dritschel 2012, Srinivasan and Young 2012, Vasavada and Showman 
2005). The reason for studying a model with all these approximations is that it is more 
tractable to analysis and interpretation. With L d = oo, the equation of motion is 
called the barotropic vorticity equation, while with finite L d it is called the equivalent 
barotropic vorticity equation. 

The (equivalent) barotropic vorticity equation is given by 

d t w + v • Vw + j3d x if — f + D, (2.2) 

where 

w = V 2 V> — (2.3) 

Here w is the vorticity, is the stream function, v = z x V/3 is the horizontal velocity. 
The deformation radius L d plays the same role here as the plasma sound radius p s 
plays in the Hasegawa-Mima equation. 

2.3 Unified Equation 

The rnHME and the (equivalent) barotropic vorticity equation can be unified in a 
single equation. In the following we also take an explicit form for the forcing / and 
dissipation D. The unified quasigeostrophic-mHME is given by 

d t w + v • \7w + [dOxf) — f + D, (2.4) 


where 


w = (2.5) 

V 2 = V 2 - aLf. (2.6) 

The nonlinear advection term v ■ Vw is sometimes written as the Poisson bracket 
where 


{A, B} = (d x A)(d y B) - (d y A)(d x B) = z x VH ■ VB. (2.7) 

We take the external forcing / to be white noise forcing £ and the dissipation operator 
to be 

D = -fiw - i/(-l) h V 2 \u, (2.8) 
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where /j is the scale-independent friction and v is the viscosity with hyperviscosity 
factor h. To recover the QG barotropic vorticity equation, set a — 1. To recover 
the physics of the mHME, set a = 1 for DW modes (k x 7 ^ 0) and &zf = 0 for ZF 
modes [k x = 0), and set L d = p s = 1. To get back to the plasma physics notational 
conventions, make the substitutions 


(x,y,'ip 1 /3,v x ,v y ) ^ {-y,x,<j>,K,-v y ,v x ). (2.9) 

In the plasma context, /3 represents the density gradient and must not be confused 
with the ratio of plasma pressure to magnetic pressure. 

In the interest of full generality, the rest of this thesis will make use of the unified 
framework. As discussed in Chapter 1, the framework sticks to the conventions of 
the geophysical literature. Unless specified otherwise, all figures will be presented 
with just the geophysical parameters in the limit L d = 00 rather than the plasma 
parameters (L d — p s — 1, olzf — 0). This is done for simplicity as there are many 
qualitative similarities in the ZF behavior in the two cases. 

2.3.1 Symmetries 

Neglecting the random forcing for a moment, we examine the relevant symmetries of 
(2.4). We assume that any symmetries of the equation are not spoiled by boundary 
conditions (e.g., take an infinite system or one with periodic boundary conditions). 
These symmetries are 

1 . x —> x + a (translational symmetry in x) 

2 . y -7 y + a (translational symmetry in y) 

3 . {y, w, ip} — )■ {—y, —w, —y} (reflection symmetry in y ) 

In other words, if w(x,y,t) is a solution, then the symmetries give us other solutions: 

1 . w(x, y , t ) = w(x + a, y, t ) 

2 . w(x, y, t ) = w(x, y + a, t) 

3 . w(x, y, t) = ~w(x, -y,t ) 

Note that since v x = dy'ip and v y = we see that the third symmetry implies 

that the solution and its symmetric partner have the same value of v x , and hence, any 
jets take the same form. In other words, on the /3-plane, eastward and westward are 
fundamentally distinguishable, whereas northward and southward in some sense are 
equivalent. There is no requirement for jet motion to be symmetric in the eastward 
and westward direction. The east-west symmetry is broken by the planetary rotation. 
In the context of plasma, the analogous statement is that there is a physical difference 
for flow in the ion and electron diamagnetic directions. 

When the random forcing is taken into consideration, exact translational sym¬ 
metry is spoiled (as is the reflection symmetry). However, we assume the forcing is 
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statistically homogeneous in space. Mathematically, this means that its statistics are 
independent of position. Then, we can still say that (2.4) satisfies the symmetries 
listed above statistically. Naively, one might expect the turbulence that results from 
a solution to (2.4) to be statistically homogeneous. This is not always the case, as we 
shall see. 

2.3.2 Nonlinearly Conserved (Quadratic) Quantities 

The 2D equation (2.4), like the 2D Navier-Stokes equation, possesses two quadratic 
quantities which are conserved by nonlinear interactions. These are the energy and 
enstrophy. The average energy density is given by 

£„ = tV / dxdy^-HV^f + iPLJ 2 }, (2.10) 

x-^y J 

where the two terms account for kinetic and potential energy. This can be rewritten 
in another form as 

If 1 

E a = -J^Y~ / dxdy-wip. (2.11) 

x y J 

The average enstrophy density is given by 

W a = —^— [ dxdy^-w 2 (2.12) 

J- J x-L J y J ^ 

2.3.3 Fundamental Dimensionless Parameter 

A fundamental dimensionless parameter controlling the zonal flow dynamics is 

(Danilov and Gurarie 2004) 

7 = £ 1 / 4 ^1/V 5 / 4. (2.13) 

This parameter is related to the zonostrophy index Rp by 7 ~ (Galperin et al. 
2010). Also, if one were to modify the definition of the small parameter a defined 
by Bouchet et al. (2013) such that the normalization length scale is the Rhines scale 
Lr rather than the size of the domain, then 7 = a -1 . Since it has been shown by 
Bouchet et al. (2013) that it is a (and ct 1//2 ) that naturally appear in the normalized 
equations of motion, we opt to use 7 instead of Rp as the descriptive parameter. 
When Ld = 00 , 7 is essentially the only independent dimensionless parameter in the 
problem. 

2.4 Quasilinear Approximation 

2.4.1 Definition 

We restrict ourselves to the quasilinear (QL) approximation of this system. Let us 
be very precise about what we mean by the QL approximation. Given some kind 
of averaging procedure, one can decompose a field into a mean plus a fluctuation. 
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Then one can write down separate equations of motion for the mean and the fluc¬ 
tuation. By the QL approximation, we mean that within the fluctuation equation 
of motion, the fluctuation self-nonlinearities (nonlinear terms involving only the fluc¬ 
tuation) are neglected. Interactions between the fluctuation and the mean field are 
retained. The equation of motion for the mean is not approximated. This definition 
is consistent with classical usage in plasma physics ( irummond and Pines 1962, Ivrall 
and Trivelpiece 1973, Vedenov et al. 1962). 

At the moment, we define our average to be a zonal average. The zonal mean of 
a quantity w is given by 


_ 1 f Lx 

w(y) = — / dxw(x,y). (2-14) 

L x Jo 

The zonal average is the conceptually simplest route, requiring the fewest number 
of assumptions, to the desired result. A substantial discussion of different types of 
averages will be given in Section 6.2. The fluctuation, or deviation from the zonal 
mean, is given by w' = w — w, and is referred to as an eddy quantity. We make this 
semantic distinction because a zonal mean quantity is likely to not fluctuate much if 
many independent correlation lengths of the turbulence have been averaged over. We 
assume the eddy quantities contain the turbulent behavior. 

To illustrate the QL approximation, we temporarily set forcing / and dissipation 
D to zero. We decompose the flow field into a zonally symmetric part (the zonal flow) 
and the residual (the eddies or turbulence). Equation (2.4) can be decomposed as 


d t w + v' ■ Vw' = 0, (2.15a) 

d t w' + v ■ Vw' + v' ■ Vw + v' ■ Vw' — v' • Vw' + (3d x ip' = 0. (2.15b) 

No approximation has been made thus far. The QL approximation involves neglecting 
the eddy-eddy nonlinearity within the eddy equation. The QL system is 

d t w + v' • Vw' = 0, (2.16a) 

d t w' + v ■ Vw' + v' ■ Vw + (3d x i\)' = 0. (2.16b) 

More explicitly, and with forcing and dissipation restored, the QL system is 

d t w' + {UV 2 + $ - [( d 2 y - L~ 2 )U}}d x ^ = £ - fiw’ - u(-l) h V 2h w ', (2.17a) 

[9t + // + o(-l) h d 2h ] (1 - a ZF L- 2 d~ 2 )U{y ) + d v v^ = 0, (2.17b) 

where 

U(y) = -d$ (2.18) 


is the zonal-mean zonal velocity. We have assumed that the zonal mean of the forcing 
is zero. If desired, one could easily allow for different dissipation rates on the zonal- 
mean quantities. 

The QL approximation does not affect the conservation of the quadratic quantities 
by the nonlinear interactions. This can be easily seen from the Fourier-space point 
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of view. Each triad interaction individually conserves these quantities, and the QL 
approximation amounts to removing some of these triad interactions. 

One issue to keep in mind is that the QL approximation breaks the material 
conservation of potential vorticity (PV). Potential vorticity, a scalar held defined by 
q = p _ 1 u? a • V6 1 , where p is the fluid density, u a is the absolute vorticity, and 6 is 
the potential temperature, is a critical quantity (Pedlosky 1987). Many quantities of 
interest can be derived from the PV, a concept known as the PV invertibility principle 
(McIntyre 2008). Furthermore, q is conserved following the how. Conservation of PV 
relies on the combination of the eddy-eddy interactions and the eddy-mean interac¬ 
tions, so the neglect of eddy-eddy interactions in the QL approximation breaks PV 
conservation. As a result, the QL system may lose certain physics that are based on 
the conservation of PV ( Dritschel and McIntyre 2008). 

2.4.2 Motivation for Using the QL Approximation 

Srinivasan and Young (2012) have shown that the QL system exhibits many of the 
same basic zonal jet features as the full nonlinear (NL) system, including the formation 
of stable jets and merging jets. With periodic boundary conditions, the equation of 
motion enjoys translational symmetry in both the x and y directions. As a parameter 
is varied, the simulations suggest a spontaneous breaking of statistical homogeneity 
in the y direction. At large /i (small 7 ), the NL system in Figure 2.1(a) and the QL 
system in Figure 2.1(d) do not exhibit steady ZFs, so the behavior is statistically 
homogeneous. At small // (large 7 ), both the NL system in Figure 2.1(b) and the 
QL system in Figure 2.1(e) do exhibit steady ZFs, implying a breaking of statistical 
homogeneity in the y direction. We also observe that, in both the NL and QL systems, 
simulations that differ only in initial conditions and realizations of the random forcing 
can display different numbers of jets [Figure 2.1 (b,c,e,f)]. I 11 addition to these features, 
both NL and QL exhibit merging jets, evident in Figure 2.1 (c,f). 

Our motivation in adopting the QL approximation is not because we believe it to 
be quantitatively correct, but rather because the QL system apparently retains the 
necessary ingredients that lead to the rich behavior of ZF formation. The QL system 
may provide insight into the more realistic models, and the advantage, of course, is 
that the QL system is far more tractable. The phenomena described above will all 
be explained analytically within the QL approximation. 

Separate from our motivations for using the QL approximation, Bouchet et al. 
have argued that in the regime of large 7 the flow becomes predominantly zonal and 
the QL approximation becomes rigorously valid ( louchet et al. 2013). Our present 
study examines the regime in parameter space where 7 is not asymptotically large, 
for it is in this regime where ZFs are born at low amplitudes from turbulence. 

The QL approximation has also been used by Herring in the study of thermal 
convection ( [erring 1963), where the only nonlinear interaction retained was between 
a horizontally-averaged temperature and the fluctuating temperature and velocity; 
nonlinear interactions between the fluctuating quantities were discarded. At large 
Rayleigh number, this approximation was able to reproduce some of the qualitative 
features observed in experiments. 
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Figure 2.1: Space-time diagrams of zonal flow in DNS of QG. Top: NL simulations at 
(a) n = 0.08 (no steady jets), (b) /i = 0.02 (8 jets), and (c) // = 0.02 (7 jets). Bottom: 
QL simulations at (d) /j = 0.29 (no steady jets), (e) /i = 0.08 (7 jets), and (f) /i = 0.08 
(6 jets). The only differences between (b) and (c) and between (e) and (f) are the 
choice of initial conditions and the realization of the random forcing. Merging jets can 
be seen in (c) at t ~ 200 and in (f) at t « 30. Our numerical simulations described here 
and elsewhere are pseudospectral, typically using a resolution of 256 x 256 (Orszag 
1969, Trefethen 2000). We use ETDRK4 as our timestepping algorithm ( box and 
Matthews 2002, Kassam and Trefethen 2005). We dealias using the 2/3 rule (Boyd 
2001, Orszag 1971). 

2.5 CE2 

The eddy quantity w' fluctuates rapidly in space and time. Averaging over these 
turbulent fluctuations enables one to work with smoothly varying functions. Such 
statistical approaches provide one path to gaining physical insight. Sometimes sta¬ 
tistical turbulence theories strive for quantitative accuracy, which requires rather 
complicated methods ( vrommes 2002), but we eschew those methods here because 
they are not required for investigation of the QL system. 

We consider an average of the QL system (2.17). The resultant framework is called 
CE2, or the second-order cumulant expansion. (If one performs a cumulant expansion 
of the original equations and truncates all cumulants higher than second order, one 
reaches the same equations.) Derivations can be found in Farrell and Ioannou ( 1003), 
Marston et al. (2008) though we follow Srinivasan and Young (2012) because there are 
advantages to that formulation. The full derivation can be found in Appendix A, but 
we give here a brief overview of the procedure. One defines the two-point, one-time 
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correlation function of vorticity using a zonal average as 


1 f Lx _ 

W(x,y 1 ,y 2 ,t) = — dx\ x w'(xi,y 1 ,t)w'(x2,y2,t), (2.19) 

Jo 

where L x is some averaging length, the integration is over the sum coordinate x = 
\{x\ + x 2 ), and the difference coordinate x = X\ — X 2 is held fixed. The correlation 
function T of stream function can be defined similarly as 

1 f Lx 

^(x,yi,y 2 ,t) = — / dx\ x ip , {x 1 ,yi,t)'ip , (x 2 ,y 2 ,t). ( 2 . 20 ) 

Jo 

One finds an evolution equation for W by taking a time derivative of (2.20), substitut¬ 
ing the expression for u/ from (2.17a), and performing the average. Under an ergodic 
assumption, the zonal average is equivalent to a statistical ensemble average, and the 
stochastic forcing can be averaged to a deterministic quantity. Then one performs a 
linear coordinate transform to the sum and difference variables y — \{y\ + 2 / 2 ) and 
y = yi — V 2 - In the ZF equation (2.17b), the Reynolds stress term can be related 
to T. The final equations are 1 

d t W(x, y\y,t) + (U + - U.)d x W - (u" + - V"_) (V + l -d*\ 


- [2/1 - (u" + + u"_)]dyd y d x ^ = F(x, y) - 2 yW - 2 uD h W, (2.21a) 

[d t + n + u(~l) h d^ h ]7U(y,t) + dyd y d x ^( 0,0 | y,t) = 0 , ( 2 . 21 b) 

where U(y, t ) is the ZF velocity, and 

U± = U(y±±y,t), (2.22) 

U" ± = U'l- & ZF LfU ± , (2.23) 

V 2 = V 2 - Lf = dl + dl - Lf, (2.24) 

7=1- a ZF Lfd= 2 , (2.25) 


F(x, y ) is the covariance of the external forcing, and Dh is the hyperviscosity operator, 
given by 


Dh = ( ^ 1) ^ 






(2.26) 


In (2.21b), the notation d y d x ^(0, 0 | y,t) implies that the partial x and y derivatives 
are taken first, and then the result is evaluated at x — y = 0. It can be shown from 

1 To transform to the conventional plasma coordinates and notation, it follows from (2.9) that 

one needs to make the substitution (a;, y , y, /3, U) 1 —>• (—y, x, x, k , —U). 
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the definitions that W and \k are related by 

W(x,y\y,t)= + d y dy + ^d (V - d y dy + T(x, y \ y, t). (2.27) 

The use of the sum and difference coordinates x, y, y allows the structure of the 
theory and especially of the bifurcation to be more easily understood than in the 
original coordinates. In the new coordinates x and y represent two-point separations 
and y represents the two-point average position. If the turbulence were homogeneous, 
there would be no dependence on y. 

The only assumption necessary for CE2 to be an exact description of the QL model 
is ergodicity in the zonal (x) direction, such that a zonal average is equivalent to an 
ensemble average. No other assumptions are required because the QL model neglects 
the nonlinear eddy-eddy term that would give rise to a closure problem. Alternatively, 
instead of the QL-based derivation, CE2 can be regarded as a truncated statistical 
closure of the NL model (Farrell and Ioannou 2003, 2007, Marston et al. 2008, Tobias 
et al. 2011, Tobias and Marston 2013). However, we prefer the former interpretation. 

CE2, like the QL system, exhibits merging jets (Farrell and Ioannou 2007). Since 
CE2 is deterministic, if the system approaches a stable steady state then merging and 
branching of jets can only occur transiently. Once the stable equilibrium is reached, 
the system is stuck there and no more dynamical behavior can occur. However, if 
the QL system is not fully ergodic, then CE2 is not an exact description of it and 
dynamical behavior like merging or branching can persist even in a statistically steady 
state (Bouchet et al. 2013, Farrell and Ioannou 2003). Though ergodicity is often 
a useful idealization, lack of complete ergodicity is to be expected in any physical 
system. 

Historically, CE2 was first studied by Farrell and Ioannou (2003) under the name 
Stochastic Structural Stability Theory, or SSST. Independently, Marston et al ( 1008), 
Tobias et al. (2011) described the second-order cumulant expansion and called it 
CE2. Later, Srinivasan and Young (201 ) also independently derived CE2 from the 
quasilinear approximation and pointed out that SSST and CE2 are mathematically 
identical. They opted to use the CE2 label, and we stick with the CE2 name for 
continuity. Recently, the acronym for Stochastic Structural Stability Theory was 
rebranded from SSST to S3T (Constantinou et al. 2013). 

2.5.1 Symmetries of the CE2 Equation 

The CE2 equations (2.21) inherit important symmetries of translation and reflection 
from the symmetries of the dynamical equation (2.4). First, we note that because 
of the form of the expression 2/3 — \U (y + \y) + U (y — |x/)], there can be no 
symmetry that changes the sign of U. Second, we examine how the two expressions 
U(y + \y) “ U(jj — \y) and U(jy + \y) + U(y — \y) behave under reflections. (The 
same expressions with U" instead of U behave in the same manner.) 


If we change y —> —y and let U(y) = U(—y), then for the first expression, 


U(y + \y) - U(y » \y) = U(-y - \y) - U(-y + \y) (2.28) 

= -U(jj+ \y) +U(jj- \y) (2.29) 


where in the second line we have let y = —y. Here, we note that the transformation 
has induced a minus sign. On the other hand, for the second expression, 

U(y + \y) + U(y - \y) = U{-y -\y)+U{-y + \y) (2.30) 

= U(v+\y)+U§-\y) (2.31) 


the transformation has not induced a minus sign. 

If we change y —» — y, y —> —y, still with U(y) = U(—y), then for the first 
expression, 


u(y + \y) ~ U{y ~ \y/ 2) =u(-y- \y) - U(-y + \y) (2.32) 

= u{v+\y)-U§-\y) (2.33) 

where in the second line we have let y = —y. Here, the transformation has not 
induced a minus sign. Similarly, for the second expression, 

U(y+b)+U(y- | y) = U{-y - \y) + U(-y + I y) (2.34) 

= U(jj+ \y) + U(jj- \y) (2.35) 


the transformation has not induced a minus sign. 

With the above relations, we can see that the equations have the following sym¬ 
metries: 


y->y + 6y, (2.36a) 

x, y —» —x, —y, assuming F(x,y ) = F(—x,y), (2.36b) 

y, y-> -y, -y, assuming F(x,y) = F(x, -y), (2.36c) 

x, y —> — x, — y. (2.36d) 


In other words, if {W(x,y \ y,t),U(y,t)} is a solution, then the symmetries give us 
other solutions: 


{W(x,y\y + Sy, t),U(y + Sy, t)}, (2.37a) 

{W(-x , y | -y,t),U(~y, t)}, (2.37b) 

{W(x, -y | -y, t), U (-//, t)}, (2.37c) 

where 5y is some constant translational shift. The symmetry (2.36d), dubbed the 
exchange symmetry, does not give a new solution because it is always obeyed by 
the correlation function such that W(x,y \ y,t ) = W(—x,—y \ y,t) ( 

Young 2012). Equation (2.36) gives the symmetries obeyed by the equations. If all 
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of the reflection symmetries are obeyed by the solutions, then one has that U is even 

in y, 

U(y,t) = U{—y,t), (2.38) 

and three relations for W: 

W(x,y\y,t) = W(-x, -y\y,t) = W(x, -y \ -y, t) = W{-x, y \ -y, t). (2.39) 

These are the symmetries in real space. We can also state what the corresponding 
symmetries are in Fourier space. It is not difficult to see that a reflection symmetry in 
real space corresponds to a reflection symmetry in Fourier space, which comes directly 
from the definition of the Fourier transform. Suppose there is some equation in x, 
and that f(x) and f(x) are both solutions, where f(x) = f(—x). Let f(k) = iF[f(x)] 

and J(k) = J r [/( x)]. Then both f(k) and f{k) will be solutions in Fourier space, 
and they will be related by 


fik) = T[f(x)\ 

-h'"* 7ix) 

= 

Thus, the reflection symmetries correspond to 


if y is kept in real space, or 


k x ,y-> -k x , -y, 

(2.40a) 

k y ,y^ —k y , —y, 

(2.40b) 

kx, ky y k x , ky, 

(2.40c) 



(2.41a) 

ky, fcy y ky , ky-; 

(2.41b) 

kxi ky y k Xl ky-, 

(2.41c) 


if y is also transformed to Fourier space. 

We also point out that as a result of the exchange symmetry, the mixed real 
spaceTFourier space quantity W(k | y, t ) is purely real and must satisfy W(k \ y,t) = 
W(-k\y,t). 
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2.5.2 Nonlinearly Conserved Quantities 

The average energy density and enstrophy density are conserved by nonlinear inter¬ 
actions in the QG-mHME and its quasilinear variant. Accordingly, they are also 
conserved in CE2. We derive here the formulas for the energy and enstrophy density 
within the CE2 description. 

First we split the average energy density, given in (2.11), into a contribution from 
zonal and eddy contributions: 


11 f y _ - 11 

Ea = ~ 2 L~ ' dyw (y^^ ~ n 


J y Jo 

= E a . Z F + -^ajeddy- 


2 L X Ly J Q 


dx / dyw'(x,y)^'(x,y ) 


(2.42) 


For the ZF contribution, recall that U(y) = — dy ip and w = — dylU , so that 


1 1 


Ea-zF = -^-j~l dyw(y)^(y) 


y Jo 


1 1 


= 2 Y' dy ^ dyIu ^ 


J y J o 


ll r h v 

= n~r / dyU(y)IU(y). 


2 L 


y Jo 


(2.43) 


For the eddy contribution, we first define a symmetrized quantity 


1 1 


E s (x,y \ y) = j dx\ x [w\x 1 ,y 1 )ilj , (x 2 ,y 2 )+ i/j , (x 1 ,y 1 )w'(x 2 ,y 2 )]. (2.44) 


X JO 


Using the same techniques as in Appendix A, we find 


E s (x,y | y) 


+ i d v *,y I y)- 


(2.45) 


The energy density E a is obtained by setting x\ = x 2 and y\ = y 2 , i.e., taking x = 0 
and y = 0, then integrating over y: 

E a - e ddy = f dyE s ( 0,0 | y). (2.46) 

1 L y JO 

Similarly, the enstrophy density can be split into ZF and eddy contributions, with 


W a = W a , ZF + W a ;eddy 5 (2.47) 

W a;ZF = ~ [ " dy [dyIU(y)} 2 , (2.48) 

z L y Jo 

fU^eddy = \ [ dy PU(0,0 | y). (2.49) 

1 L y Jo 
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2.5.3 Wigner-Moyal Formalism 

The Wigner-Moyal formalism, which has been used in studies of wave physics in 
inhomogeneous media ( lall et al. 2002), is basically equivalent to CE2. The Wigner 
distribution function, assuming an appropriate average is used in its definition, is 
closely related to the CE2 correlation function W: they are both the two-point, one¬ 
time, second-order correlation of fluctuations. The Wigner-Moyal equation, which 
describes the evolution of the distribution function, is the analog of the CE2 equation 
(2.21a). The Wigner-Moyal formalism has also been used as the starting point for a 
few calculations involving zonal flows (Mendonga et al. 2014, Mendonga and Benkadda 
2012, Mendonga and Hizanidis 201 ). Those papers perform some analyses similar to 
what is in this thesis, but they made several further approximations without stating 
regimes of validity. In contrast, this thesis provides a deep understanding of the 
theory without further approximation and also frequently compares analytic results 
with numerical results to ensure correct understanding. 

2.5.4 Wave Kinetic Equation 

Other previous studies of zonal flows have used a wave-kinetic framework of inhomo¬ 
geneous turbulence (Diamond et al. 2005, Dyachenko et al. 1992, Krommes and Kim 
2000, Krommes and Parker, Manin and Nazarenko 1994, Smolyakov et al. 2000b). 

The wave-kinetic formalism, like CE2, describes fluctuations using a second-order, 
two-point correlation function. In these studies, the wave-kinetic formalism is re¬ 
stricted such that the length scale of the inhomogeneity must be much longer than 
the small scales of the turbulence. A traditional viewpoint is that the wave kinetic 
equation describes turbulence as wavepackets that propagate through an inhomoge¬ 
neous medium. CE2, on the other hand, is an exact description of the QL equations 
and makes no approximation or restriction on length scales. 

The disparate-scale asymptotic limit of CE2 recovers the wave kinetic equation. In 
(2.21a), assume d y <C d y , Taylor expand the terms U(y ± \y), and Fourier transform 

(x,y) —* (k x ,k y ). After switching to using A^(k | y) — (1 — azpk Ly 2 )W( k | y) as 
the dependent variable rather than W{ k | y ), the disparate-scale form of CE2 takes 
on wave-kinetic form. 
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Chapter 3 

Zonostrophic Instability and 
Beyond 


This chapter develops the physics at the core of this thesis. First, we review zonos¬ 
trophic instability (ZI). In ZI, a statistically homogeneous turbulent state is unstable 
to coherent, zonally-symmetric perturbations. These perturbations grow into zonal 
flows (ZFs). Since ZI is an instability of a turbulent, time-dependent but statistically- 
steady state, analysis of it requires a statistical formalism. CE2 provides the sim¬ 
plest such formalism, and indeed, the instability was discovered through the CE2 
framework. The fundamental dynamical equations are too complicated for analytic 
progress. 

ZI has been explored numerically and calculated analytically in detail. We provide 
a thorough review of the analytic calculation and explore certain limits of the dis¬ 
persion relation. Our calculation mildly generalizes previous work because we allow 
for finite deformation radius (or Larmor radius) Lj as opposed to infinite L^ and we 
allow for viscosity in addition to a scale-independent drag. 1 

We also draw a connection between ZI and generalized modulational instability. 
By modulational instability we mean the instability of a primary mode to a secondary 
mode. We find that with the ZI dispersion relation from CE2, we can recover as a 
special case a previously-derived dispersion relation of modulational instability. This 
discovery suggests that CE2 may be useful for future investigations of modulational 
or secondary instabilities or generalizations thereof. 

We then extend analytic understanding of ZI beyond a linear stability calculation 
into the regime of nonlinearly interacting zonal flows and turbulence. We connect 
the generation of zonal flows to the large body of literature of pattern formation. 
At a basic level, zonal flows appear in a spontaneous symmetry-breaking bifurcation 
where the broken symmetry is statistical homogeneity in space. The mechanism of 
the symmetry breaking is ZI. 

We perform a bifurcation analysis, which yields numerous insights. We construct 
explicit solutions to the nonlinear CE2 equations, and we discover important and 

1 Srinivasan and Young (2012) have some calculations in an appendix that include viscosity, but 
there is an error in the way the wavevector dependence of viscosity terms is treated. 
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unexpected features. First, we find that the zonal flow wavelength is not unique. 
Many wavelengths allow a steady-state solution to the equations. Second, only some 
of these wavelengths correspond to solutions which are stable. Unstable wavelengths 
must evolve to reach a stable wavelength; this process manifests as merging jets. 
Consequently, we are able to provide a theoretical basis for the well-known merging 
of jets along with a simple PDE that demonstrates the behavior. Furthermore our 
work links the merging of jets to the large body of research of defects, providing new 
avenues for research into jet dynamics. Our results provide a substantial theoretical 
foundation for further understanding of turbulence and zonal flows. 

This chapter is structured as follows. Section 3.1 introduces a phenomenological 
model of the bifurcation. As a zero-dimensional model, it orients the reader before 
the plunge into the full problem with its complexity of spatial dependence. A review 
of ZI is provided in Section 3.2. This calculation shows that a state of turbulence 
without zonal flows can be unstable to zonal flow perturbations. We discuss how ZI 
relates to modulational instability in Section 3.3. Then, in Section 3.4 (with details 
in Appendix C), we perform a full bifurcation analysis into the regime of nonlinearly 
interacting eddies and zonal flows. 

3.1 Phenomenological Bifurcation Model of Zonos- 
trophic Instability 

A zero-dimensional phenomenological model illustrates some of the key features of ZI 
and the bifurcation to a state with ZFs (Parker and Krommes 201 1). The system is a 
variant of another treatment which models the appearance of shear flows in the L-H 
transition in plasmas ( Diamond et al. 1994). However, the model we present more 
closely mirrors the structure and behavior of the CE2 equations. The model includes 
three interacting degrees of freedom: the homogeneous, or spatial average, part of the 
fluctuation covariance IIpp the inhomogeneous, or deviation from the spatial average, 
part of the fluctuation covariance W % ; and the ZF amplitude (not covariance) z. Both 
Wi and z may be positive or negative. The model is given by 


Wh = —^Wh - aWiZ + F, (3.1a) 

W i = -ji.W i + riW h z, (3.1b) 

z = —vz + aW t . (3.1c) 

The structure of the model reflects that of the CE2 equations (2.21) in several 
ways. To affect the homogeneous part of the turbulence, the ZF interacts only with 

the inhomogeneous part. Similarly, the ZF interacts with the homogeneous part 

to affect the inhomogeneous part. (CE2 also contains an interaction between the 
ZFs and the inhomogeneous part to affect the inhomogeneous part; this is neglected 
here, as in ZI analysis.) Finally, it is the inhomogeneous part of the turbulence 
that is responsible for driving steady ZFs. The above model neglects the eddy self- 
nonlinearities as in CE2. The appearance of the same coefficient a in Wh and in 
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z reflects the conservation by the nonlinear interactions of an energy-like quantity 
Wfr + \z 2 . We take all of the coefficients /i, u, a , ?/, F to be positive. 

The model allows a homogeneous equilibrium, in which forcing F is balanced 
by dissipation /i, and for which W t and z are zero. The homogeneous equilibrium 
is unstable if Fr\ajp 2 v > 1. Increasing the forcing or decreasing the dissipation 
tends to make the homogeneous equilibrium more zonostrophically unstable, which 
is characteristic of the more rigorous analysis. 

When the homogeneous equilibrium goes unstable, it connects to an inhomoge¬ 
neous equilibrium at Wh = pu/rja, W? = ( u/a 2 )(F — p 2 u/r}a), z = aWijv. This 
new equilibrium is stable (when it exists), which can be seen by constructing the 
eigenvalues graphically from the characteristic polynomial. Furthermore, there are 
actually two symmetric solutions, with either sign of £ and W % . There is thus a su¬ 
percritical pitchfork bifurcation; this feature is also present in the complete model, 
but the discrete z —>■ — z symmetry becomes a continuous symmetry associated with 
translational invariance. 

The model demonstrates some of the qualitative features of ZI, although in sim¬ 
plifying it we have tossed out spatial dependence. Spatial dependence makes the 
problem both immensely more complicated and immensely more interesting. The 
CE2 equations contain the full spatial dependence. Detailed analysis of ZI and be¬ 
yond proceeds in the next few sections. 

3.2 Zonostrophic Instability 

In this section we review ZI, for which substantial understanding has been recently 
obtained (Bakas and Ioannou 2011, Srinivasan and Young 2012). To give a brief 
overview, ZI refers to an instability where a state of homogeneous turbulence without 
ZFs can be unstable to ZF perturbations. In the regime where ZI is present, inho¬ 
mogeneous turbulence results. The instability as well as the nonlinear growth and 
saturation can be handled self-consistently within the CE2 framework. This section 
is devoted to the study of the instability of the homogeneous equilibrium, with later 
sections handling the nonlinear saturation. 

We examine the homogeneous equilibrium of the CE2 equations, which has no 
zonal flows. We calculate the linear response of the equilibrium to zonal perturbation. 
Much analytic progress is possible, which provides substantial insight. Although the 
final dispersion relation must be solved numerically, it can be reduced to a single 
nonlinear equation. In some regimes of parameter space the equilibrium is unstable, 
and the instability has been named zonostrophic instability. This instability has been 
studied analytically in detail (Srinivasan and Young 2012) and aspects of it have also 
been examined numerically (Bakas and Ioannou 2011, Farrell and Ioannou 2007). 

As a control parameter p is varied, the homogeneous state becomes zonostroph¬ 
ically unstable (Farrell and Ioannou 2007, Srinivasan and Young 2012). Physically, 
ZI occurs when dissipation is overcome by the mutually reinforcing processes of eddy 
tilting by ZFs and production of Reynolds stress forces by tilted eddies. The in¬ 
stability eigenmode consists of perturbations spatially periodic in y with zero real 
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frequency (Srinivasan and Young 2012), so that zonostrophic instability arises as a 
Type I s instability (Cross and Hohenberg 1991) of homogeneous turbulence. 


3.2.1 CE2 Homogeneous Equilibrium 

A homogeneous, steady-state solution of the CE2 equations always exists, arising 
from a simple balance between forcing and dissipation. This solution is 

W H = (2// + 2uD h )~ 1 F, (3.2a) 

U = 0, (3.2b) 


where the H subscript denotes homogeneous. From (2.27) it is easy to relate Wh and 

h'■ _ 

W H (x, y ) = V V H (x, y )• (3.3) 

We can also give the result in Fourier space by applying the continuous Fourier 
transform 

W(k x ,k y ) = J dxdye- ik * x e~ ik y y W{x,y). (3.4) 

For the homogeneous equilibrium with ky = 0, we have 2vDh —> 2uk 2h , where k 2 = 
k 2 + k 2 . This gives for the homogeneous equilibrium 


1 ^Ii{kxi ky ) 


F ( k x , ky ) 

2 (/r + uk 2h ) 


and 

( k x , ky ) A T// (k x , ky '), 


where 

k 2 = k- + Lf. 


(3.5) 

(3.6) 

(3.7) 


3.2.2 Linearization about the Homogeneous Equilibrium 

The homogeneous equilibrium is linearly stable in a certain regime of parameters. To 
determine its stability one calculates the dispersion relation corresponding to ZI. One 
considers perturbations about the equilibrium in (3.2). The derivation given here 
closely follows that given by Srinivasan and Yonng (2012). Because the equilibrium 
is independent of y and A, the y and t dependence of the perturbations can be Fourier 
transformed. The fields are written as 

W(x, y \y,t) — W H (x, y) + 5W{x , y)e xt e i (3.8a) 
U(y,t) = dUe xt e ig y, (3.8b) 

where q is the ZF wavenumber and A is the eigenvalue. 
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We now substitute the perturbations into (2.21a) and (2.21b) and linearize. We 
use that 


U± = U(y ± y/2) 6Ue xt e iqy e ±iqy/2 , (3.9) 

U'± —q 2 U±, (3.10) 

U± ->• —q 2 U± - a ZF LfU± = -q 2 SUe xt e iqy e ±iqy/2 , (3.11) 

where q 2 = q 2 + &zfL^ 2 . The linearized equations are 

AhW + 5U(e iqy/2 - e~ iqy/2 )d x W H + q 2 5U{e iqy/2 - e~ iqy/2 )V' 2 d x ^ H 

- i2/3qd x d y 5 T = -2 y6W - 2 uD h 6W, (3.12a) 

(X+H + vq 2h ) (l + azFL~ d 2 q~ 2 )5U = —iqd x d y 5^{0 , 0). (3.12b) 


Note that we can express 1 + azFL^ 2 q~ 2 = q 2 /q 2 ■ If for wavenumber q , A is an 
eigenvalue with eigenvector (5W,SU), then for wavenumber —q, X* is an eigenvalue 
with eigenvector (5W*,5U*). 

It is convenient to Fourier transform in both x and y as well. For the perturbations 

—2 —2 

(2.27) becomes, with d y —> iq, d y —> ik y , and V —» —k , 


5W(k x ,ky) = hlh 2 _8*(k x ,ky), 
h 2 ± = k 2 x + (. k y ± q/ 2) 2 , 

hi = hi + T- 2 . 


We also use 


d x d y f(^Xj y) |x,y=o,o d x d( 


y {2n) 2 

1 f 


(2tt) 


J dk x dk y e ikxX e ikyV f(k x , k y )\ x , y= o,o 

dk x dk y kjzkyfi^kxi k y ). 


Thus the ZF equation (3.12b) becomes 


g 2 


— (A + y + z/g 2ft ) <57/ = iq 


dk x dky \ 2 


kxky 5^(k x ,k. 


y> 


(3.13) 

(3.14) 

(3.15) 


(3.16) 

(3.17) 


(3.18) 


Now we transform the DW equation. The second and third term of (3.12a) can 
be combined using (3.3) as 

SU(e iqy/2 - e ~ iqy/2 ) V 2 (V 2 + f)d x ^ H 


. Using the property that J r \e iqy ^ 2 f(y)\ = f(k y — q/2), the Fourier transform of 


e±i qy / 2 V 2 (y 2 +f)d x * H (x,y) 
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IS 


ik x [L d 2 + k 2 x + ( k y =f q/2) 2 ] [L 2 + fc 2 + (k y =f g/2) 2 - g 2 ]'F# (/;*, =f §g) 


= ik x h\(h 2 ^ - q 2 )^ H {k x ,k y =F \q). 

(3.19) 

Also, Dh transforms to 


D h = ^{h 2 + h + li 2 _ h ). 

(3.20) 

Equation (3.12a) then becomes 



A/i+/i!.(J^(A: x , ky) + ikJU hi {hi 

+ i2/3qk x k y 5'$ l = 


—2 / t -2 


- (A:*, ky ~ \q) ~ h 2 + {h\ 

- [2fi + u{hf + h 2 _ h )]h 2 + hi6V. 


q 2 )^ H {k x ,k y + \q) 
(3.21) 


Let 

= h± { ,l ± ~ Q 2 ) ^ h ( K , k y ± §g). (3.22) 

Rearranging slightly, the linearized equations about the homogeneous equilibrium are 


h 2 + h 2 _ (A + 2/i + i /(h 2/l + h 2 ! 1 )) + i2fiqk x k y 5^{k x , k y ) + ik x 5U{<$>~ H - <h+) = 0, 

(3.23a) 

k x h,^„ ( 3 . 23b ) 


q 2 


— (A + n + vq 2h )6U = iq j dk x dk y -^^6^{k x ,k y ). 


Here, (3.23a) and (3.23b) are exact equations for the eigenvectors. However, for given 
parameters, only certain values of A allow eigenvectors. Those are the eigenvalues. 
We can determine the values of A which give solutions by using (3.23a) to solve for 
S^/{k x , k y ) in term of SU, then substituting into (3.23b). A nonlinear equation results. 
Once we know the eigenvalues A, we can find the eigenvectors by taking some value 
for SU and using (3.23a) to give the 5W(k x , k y ). 


3.2.3 Dispersion Relation 

We now obtain the dispersion relation. First we solve for hT in terms of 5U: 


t 


— (A + n + vq 2h )8U = iq J dk x dk. 


k x k y 


-ik x {<b- H - <$> + H )SU 


(2tt) 2 h + h_ (A + 2 n + v{h 2 + + h 2 ! 1 )) + 2i/3qk x k y 

(3.24) 

We substitute this back into the equation for SU, and also rewrite in terms of 
Wh • This yields the dispersion relation 


t 


^2 (A + 11 + vq 2h ) = qA- ~ qA+, 


(3.25) 
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where 


A + = 


dk x dk y 


k x k y (l - q 2 /h 2 ± )W H (k x ,k y ± \q) 


( 2 T> 2 [A + 2 fi + u(hf + h 2 _ h )] h 2 + h 2 _ + 2 i/3qk x k y ’ 


(3.26) 


^ _ 2 

and h± = kl + (k y ± and h ± = h\ + L^ 2 . Some algebraic manipulation shows 
that A + = —A_ ( ). This is done by first noting that 

W(k x , k y ) = W(—k x ,—k y ) for any correlation function, and then letting k x —> —k x 
and k y —>• —k y in the integral for A + . 

This dispersion relation was also obtained by darnevale and Martin (1982), in a 
form that allowed for arbitrary inhomogeneities rather than only zonally symmetric 
ones. That paper did not, however, remark on the connection to the generation of 
zonal flows. 

Equation (3.25) is the general dispersion relation. Following hinivasan and Young 
(2012), we also provide the specialized results for an isotropic turbulent background 
spectrum. Although a purely isotropic spectrum is unlikely to obtain in practice when 
the beta effect is present, such an investigation helps to get a simplified dispersion 
relation, gain intuitive understanding, and isolate the physical consequences of various 
effects. 

In A_, make the transformation k' y = k y — | q . After working through the trans¬ 
formation, then dropping the prime on k' yl we End 


2gA_ = 


dk x dky 


2 qkl(k y + \q) (l - q 2 /k 2 )W H (k x , k 


( 2 T> 2 [A + 2 fi + u(hf + + k 2h )] h 2 ++ k 2 + 2ifdqk x {k y + q/2) ’ 


(3.27) 


where K 2 ++ = k 2 + (k y + q) 2 and h~ ++ = h 2 ++ + L d 2 . Also note that one can write 
h 2 ++ — k 2 = 2 q{k y + \q). Now rewrite the integral using polar coordinates, with 


—keosej), and note 


2 q ( k y + \q) = k 2 (n 2 — 2 n cos <p) , 

(3.28) 

h 2 ++ = k 2 (l — 2 n cos </> + n 2 ) , 

(3.29) 

k 2 = k 2 ( 1 + m), 

(3.30) 

h ++ = k 2 (l — 2 n cos p + n 2 + m), 

(3.31) 


where n = q/k and m = (/cL^) -2 . Assuming the equilibrium is isotropic, 
W H (k x ,k y ) = Wn(k), then after some manipulation the dispersion relation can 
be put into the form 




/(A + 2^^ , h 


fa 


Pq 1 k 


(3.32) 
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where 


S(x,V, n ,m,h) = J ^K, (3.33) 

_ (n — 2 cos 0) sin 2 </> 

jfV = -. 

{x + ??[1 + (1 — 2n COS0 + n 2 )^]}^! — 2n cos</> + n 2 + m) + i(n — 2 cos 0) sin 0 

(3.34) 

We now specialize to thin-ring forcing, where the wavevectors excited by the ex¬ 
ternal forcing are confined to a thin ring in k-space. We take 


F(k) = 47 TekfS(k — kf), 


(3.35) 


where £ is, in the case of L d —> oo, the total energy (density) input. Then, from (3.5), 
the homogeneous equilibrium is 


W H (k) 


2nekf 
H + ukj h 


S{k-kf). 


(3.36) 


Substituting this into (3.32), the integral over the delta function is trivial and we 
obtain 


^(\+Li+vq 2h ) 


£ kj(l - q 2 /k~ f ) /(A + 2/j,)k 2 f vkfk) 

P li + vkf \ J3q ’ Pq 


q_ 


0 k f L d )" 2 , 



(3.37) 


_2 

where k^ = k 2 + L^ 2 - This nonlinear equation for A involves only one integral—the 
polar integral in S —that must be computed numerically. 


3.2.4 Behavior of the Dispersion Relation 

We begin by showing some examples of the dispersion relation solved numerically. In 
each case we use the thin-ring forcing just described. We do not attempt to draw 
any definitive conclusions from the few examples we show here, but rather use them 
to get a general sense of how the dispersion relation behaves. Then, we analytically 
explore a few limits of the dispersion relation. 

Numerical Results 

First, in Figure 3.1 we show the behavior of the LHS and RHS of (3.25) (more 
precisely, the isotropic version in (3.37)) as a function of A. In the example shown, 
there is an intersection at positive A, so instability occurs. As reported by Srinivasan 
and Young (2012), numerical results indicate that all the unstable A’s are pure real, 
and we have found the same. 

Next, we plot A (q) in Figure 3.2. As a parameter such as /i is varied, the equilib¬ 
rium can go from being stable to zonal perturbations (A < 0 for all q), to having a 
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Figure 3.1: Plot of the RHS and LHS of (3.25) for real A. In this example, there is a 
positive eigenvalue solution, so zonostrophic instability occurs. 

single marginally stable mode (A = 0 at one q), to having a band of unstable modes 
(A > 0 for some q). 

Another result of interest is the neutral curve. If /i is our control parameter, then 
the neutral curve is the curve in (q,qi) space given by A (q, n) = 0. The neutral curve 
is the boundary between zonostrophically stable and unstable regions. Below the 
bottom of the neutral curve, the homogeneous state is stable. Above the bottom of the 
neutral curve, at a fixed value of /i, the homogeneous state is unstable to perturbations 
with wavenumbers q inside the neutral curve. An example of a neutral curve is shown 
in Figure 3.3 (negative ji is plotted because a neutral curve conventionally opens 
upward). 

So far we have only been concerned with real eigenvalues A. In looking to see 
whether there are any complex eigenvalues at all, including damped ones, we show 
in Figure 3.4 the residual of the dispersion relation (the difference of the RHS and 
LHS of (3.25)) as a function of A with other parameters fixed. We find that there is 
a single damped, complex eigenvalue (along with its complex conjugate). 

Analytic Limits 

We now explore the dispersion relation analytically by examining various limits. We 
explore the small q (long ZF wavelength) limit and the effect of an isotropic spectrum. 


41 












q 

Figure 3.2: Plot of the dispersion relation A (q) for several values of //. 


v=0.001 (3=1 h=4 L“ 2 =1 a =0 

r d ZF 



Figure 3.3: Neutral curve. For values of fi below the bottom of the curve, the homo¬ 
geneous state is stable. For other values of /r, the homogeneous state is unstable to 
perturbations with wavenumbers q inside the neutral curve. 
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Figure 3.4: Magnitude squared of the residual of (3.25) plotted on a logarithmic scale, 
in the eigenvalue complex plane. Along with three real eigenvalues, there is a single 
eigenvalue with nonzero imaginary part (along with its complex conjugate), which 
has negative real part. This is at q = 0.5. 



First, it is easy to find the small q limit of the general dispersion relation (3.25). 
Noting that Ii+' + h 2 ! 1 = 2 k 2h + 0(q 2 ) and h + h_ = k + 0(q 2 ), then keeping only to 
first order in q, we have 


A+ = 


dk 

W)‘ 


Co ± x q 


1 dc n 


2 dL, 


y J L 


W H (k x ,k y ) ± -q 


1 dW, 


H 


k^y 


—-2 


2 dky _ [A + 2 fi + 2uk 2h ]k 41 + 2 i/3qk x k y 

(3.38) 


where c 0 = 1 — O-zfL w Z k . Then 


gA_ — qA + = —q 2 


dk k^k y d(coWh) 

(2' 7r )" [A + 2/i + 2uk 2h ]k 4 + 2i/3qk x k y 


(3.39) 


If the dissipation terms /i and v are not small enough to be negligible, then the fd 
term in the denominator should be neglected as small in q. However, if dissipation is 
negligible, then the fd term cannot be ignored in general because A also turns out to 
be small in q. 

Second, we examine the dispersion relation for the special case of an isotropic 
spectrum. In the context of an infinite deformation radius L d , the effect of an isotropic 
background spectrum has been studied before (Bakas and Ioannou 2013b, Srinivasan 
and Young 2012). Those studies concluded that for an isotropic background, fd ^ 0 
is required for instability. Additionally, they found that for an isotropic background, 
the eddies acted on long-wavelength zonal flows as a negative hyperviscosity instead 
of negative viscosity. That is, the eddy forcing on the RHS of (3.25) behaves as g 4 
rather than q 2 at small q. In this section, we study how these results change when 
finite deformation length is allowed. 
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The dispersion relation for an isotropic spectrum is given in (3.32) and involves 
the function S(xi T], n ,m,h). For simplicity, we ignore viscosity which corresponds 
to setting r] to zero (in which case the hyperviscosity factor h drops out also). We 
explore the limit of large y, which could correspond to either small f3 or small q. 
Asymptotic expansion of 5'(y, n, m) for large y reveals interesting behavior that can 
differ for finite vs. infinite L d . 

For infinite L d (i.e., m = 0), S behaves as 2 ( inivasan and Young 2( ) 


S(y,n,0) 




n 3 

^ 8 ( 1 ^ 

1 n 2 — 1 
y 2 n 3 


n 2 ) 


+ 0( y- 5 ), 


+ 0(y" 3 ), 


n 2 < 1, 


n 2 > 1. 


(3.40) 


For small q , we recover S ~ g 4 . Additionally, we can consider the case of finite q 
but small f3. For n 2 < 1, the RHS of (3.32) goes as /3 2 , which vanishes at [3 — 0. 
Therefore, at /3 = 0 any thin ring of an isotropic spectrum with k > q has no net 
effect on the zonal flow. On the other hand, for n 2 > 1 the f3 dependence in the RHS 
of (3.32) vanishes. Thus, at f3 = 0 a thin ring with k < q has a net damping effect 
on the zonal flow. 

For finite L d , S behaves as 3 


S^y, n, m) = (4n 3 y) 1 — n 2 (—1 + m) + (1 + m) ^ — 1 — m 
+ a/[(-1 + n ) 2 + m\[( 1 + n ) 2 + m]) + O (yf 3 ) • 


(3.41) 


For small /3, the /3 dependence cancels out of the RHS of (3.32). Hence, instability 
is possible even with f3 = 0. For concreteness, one might take m — 1, for which S 
simplifies to 

S{x ' n ' 1)= ^k r 1+ M) + 0 (x " 3) • (a42) 

Additionally, the small q limit of (3.41) is 


S(x,n,m) 


n m 
X 2(1 + m ) 2 ~*” 


(3.43) 


Thus, for an isotropic spectrum and finite L d (and also n — 0), the growth rate of 
the zonal flows goes as q 2 at small q , rather than like g 4 as in the case of infinite L d . 


3.3 Connection to Modulational Instability 

Zonostrophic instability can be understood in a very general way as the instability 
of some turbulent background spectrum to a (zonally symmetric) coherent mode. 

2 Validity of this formula requires that 1 — n 2 is not too small. 

3 Validity requires that m ^ 0, because for to = 0 and n 2 < 1, the lowest order result vanishes. 
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As a special case, one can consider the background spectrum to consist of only a 
single mode. Parker and Ivrommes show that in this case the dispersion relation of 
zonostrophic instability reduces exactly to that of the 4-mode modulational instabil¬ 
ity (sometimes called parametric instability). This correspondence was first noted 
by Carnevale and Martin (1982), but they did not discuss it in the context of the 
generation of zonal flows. 

The stability of a single, primary wave p to perturbations is a problem that has 
received attention in the past (Connaughton et al. 2010, Gallagher et al. 2012, Gill 
197 1, Lorenz 1972). These calculations have used the fluctuating dynamical equa¬ 
tions such as (2.4) and not a statistically averaged system. Generally one considers 
the unforced, undamped case, for which a single wave is an exact solution of the 
nonlinear dynamical equations. Conceptually similar is the so-called secondary in¬ 
stability, where a growing, primary eigenmode gives rise to a secondary mode ( hunk 
2007, Pueschel et al. 2013, Rogers et al. 2000). If the secondary mode grows much 
faster, the primary mode is treated as a stationary background. These secondary 
instabilities are more complicated, since due to the toroidal geometry, the growing 
eigenmode has nontrivial spatial dependence. Additionally, the eigenmode is not an 
exact solution of the nonlinear equations. 

To calculate the stability of the primary wave using (2.4), in general one needs 
to retain an infinite number of coupled, perturbing modes. However, typically one 
truncates the system, for example retaining a secondary mode q and the sideband 
pair p ± q. Within this 4-mode approximation and the further assumption that the 
primary has p y = 0 such as a pure Rossby or drift wave and the secondary has q x = 0, 
the dispersion relation for 4-mode modulational instability is given by (Connaughton 
et al. 2010) 


2M 2 (1 - s 2 )(l + a 2 + /)(! + /) 2 - (s 2 + /) 

(l + /) 2 (l + s 2 + /) 2 (s 2 + /) 

where A' = pX/j3, s = q/p, f = p~ 2 L~^ , M = Vbp 3 //3, and -0 O is the amplitude of the 
background stream function. 

Some studies investigated this phenomenon by using a form of CE2 where the inho¬ 
mogeneity is assumed to vary slowly in space compared to the turbulence ( Oubrulle 
and Nazarenko 1997, Manin and Nazarenko 1994, Smolyakov et al. 2000b, Trines 
et al. 2010, Wordsworth 2009). With that assumption, the turbulence is described by 
a wave kinetic equation. The wave kinetic equation can also be recovered from CE2 
as described in Section 2.5.4. While those previous studies are limited to the regime 
of small q, the CE2 framework makes no assumption about the length scale of the 
inhomogeneity. Moreover, those previous studies did not draw a direct connection 
between the results from the statistical calculation and from the 4-mode calculation. 4 

4 One reason a connection may not have been made is that the small-g results in Manin and 
Nazarenko (1994) and Smolyakov et al. (2000b) based on the wave kinetic equation are incomplete. 
Their dissipationless (/i = 0) formulation amounts to neglecting the term 2 i/3qk x k y compared to A 
in the denominator of (3.26). But this is invalid if A ~ q 2 because the neglected term is larger than 
the retained term. For example, when specialized to a single primary mode, both papers state that 


(3.44) 
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This dispersion relation (3.44) can be recovered from CE2 and the zonostrophic 
instability dispersion relation (3.25). To precisely compare, one must carefully select 
the background spectrum Wh to correspond to a wave of stream function if; 0 . If the 
initial background amplitude of mode p is i/jq, then we write 

V>(x) = V>o (e ipx ~ iuJt + e~ ip x+iujt ) . (3.45) 

Appendix B shows that this corresponds to a one-time, two-point covariance of 
strea m fun ction 


'k H{k x , k y ) = (27t) 2 -0o [h(k - p) + 5(k + p)] . (3.46) 


From (2.27), the corresponding covariance of vorticity is given by Wn{k x , k y ) = 

—4 

k T H(k x , k y ), and thus, because of the delta functions, 

W H (k x , k y ) = (27r) 2 A[h(k - p) + <5(k + p)], (3.47) 

where we have defined A = ^q( p 2 + L^ 2 )". There are two ways of achieving this 
background spectrum. First, we could choose the external forcing to be F( k) = 
2/iWh- Since we want the dissipation term p to disappear in the final expression, 
H can be chosen to be vanishingly small, in particular smaller than the eigenvalue 
A. Alternatively, as previously mentioned we could take the external forcing and the 
dissipation to be zero, in which case any arbitrary homogeneous spectrum trivially 
satisfies the CE2 equations. This latter point of view is closer to the traditional 
stability calculations. 

Substituting (3.47) into (3.25), we find 


t 


— A = 2 qApJ 1 - =2 


r 


f 


Py + \q 


Vy-~ 2 q 


p 2 ) \^P 2 +P 2 + 2 if3qp x (p y + I q) A pip 2 + 2 if3qp x (p y - \q) 


(3.48) 

where dissipation has been neglected, = p 2 + ( p y ± q) 2 , and p\ = p\ + L^ 2 . 

When specialized to the case of a primary wave with p y = 0, the dispersion relation 
becomes 

2 A p 2 + p 2 




—A = 2 qAp x 


1 - Si 7 


(3.49) 


" ~ \ p" J 2 \ 2 p+p 4 + /3 2 q 4 p 2 ' 

Now, taking &zf = 1 to specialize to quasigeostrophic physics and introducing the 
same normalizations as used in (3.44), we obtain 


s 2 + f y = 2As 2 \'(l-s 2 )(l + s 2 + f) 

s 2 (/Vp) 2 [A' 2 (1 + s 2 + /) 2 (1 + f) 2 + s 4 ] ’ J 

for the (unmodified) Hasegawa-Mima Equation, instability occurs when p\ + L y 2 — 3 p 2 > 0, and 
that A ~ q 2 . When the /3 term is unjustifiably neglected, this result can be found from the small q 
limit of (3.48). Careful analysis shows this result also obtains in the ipo ► oo limit. But contrary to 
statements made by Connaughton et al. (2010), the wave-kinetic formalism is not restricted to that 
large-amplitude regime. If the /3 term is retained, the full answer at small q can be recovered from 
the wave-kinetic formalism. 
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Letting A 1 = p 2 A/(3 2 , after some simplification we find 


/3 _ 4 f ^A'jl — S 2 )(l + S 2 + /) — (s 2 + /) \ 

‘ V (l + f) 2 (l + s 2 + f) 2 (s 2 + f) )■ 


(3.51) 


Since A' = p 6, 0q( 1 + /) 2 //3 2 = M 2 (l + /) 2 , this exactly matches the dispersion relation 
given in (3.44). 

It may be at first surprising that the two dispersion relations agree exactly, but 
retrospectively it makes sense. The 4-wave modulational instability contains the 
primary wave p and the perturbations at wave vectors q and p ± q. From (B.5) in 
Appendix B for the correlation between the primary mode k = px and sidebands k' = 
px±gy, we see that the spatial dependence of the correlation goes as cos (px±^qy±qy). 
Upon examining the CE2 calculations, we see that the retained modes are the zonal 
flow SUe ±iqq (which corresponds to mode ±q) and the perturbations to the spectrum 
5W(k x , k y )e ±iqq . The perturbation SW(k x , k y ) is proportional to Wn{k x ,k y ± \q ), 
which is nonzero at k x = p and k y = 3n\q for the given primary mode. Therefore 
the perturbations kept within CE2 are precisely the corresponding modes kept in 
the 4-mode truncation. The CE2 instability calculation neglects higher harmonics 
of q such as e 2iqq at the linear level. These higher harmonics are precisely what is 
neglected by truncation to 4 modes instead of retaining higher sidebands. 

In the above calculation, we have shown that from CE2 we recover the 4-wave 
modulational instability in the special case of a primary wave with p y = 0 and a 
secondary wave with q x = 0. We now generalize this to show that CE2 recovers 
the 4-wave modulational instability for an arbitrary primary wave and an arbitrary 
secondary wave. 

The 4-wave modulational instability has the dispersion relation (Connaughton 
et al. 2010) 


(<? 2 + L d 2 )\ - i/3q x = ^olP x q| 2 {P 2 ~ <? 2 ) 


p 2 + — p 2 


p 2 — p 2 


(p 2 + L d 2 )(\ + iu) + i/3 ( p x - q : 


(p+ + L d 2 )(X-iu) - i/3(p x + q x ) 

(3.52) 


where p± = p ± q and oj = ~/3p x / (p 2 + L d 2 ). 

To allow for an arbitrary secondary wave within the CE2 formalism, we use the 
recent formulation of Bakas and Ioannou (2013a, c). That formulation allows for 
coherent structures of arbitrary spatial dependence rather than restricting to zonally 
symmetric structures. (The rest of this thesis is focused on zonal flows and uses the 
formulation only for zonally-symmetric structure.) Their formulation also assumed 
infinite deformation radius, though that could be modified. The dispersion relation 
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in the small forcing and small dissipation limit is (Bakas and Ioannou 2013c) 5 


A q 2 - i/3q x = 

where 

N 2{k x q y k y q x / n q x qy 


dk x dk y N 
(2tt) 2 D 


i i ^ x 

rv x “1“ ~T" 
2 


f 1 _ J// 2 ^ w H(k x , ky ), (3.53) 


% 


k v + ^ 
y 2 

+ (tfj - Qx) ( k x + y) (k y + y) k 


(3.54) 


= A/c 2 fc 2 - -~iq x !3 [k 2 + k 2 + ] + 2if3 ( k x + ^ )q x + [ky + tt )% 


(3.55) 


and k + = k + q. As before, the appropriate background spectrum to correspond with 
that of (3.52) is Wh = (27t) 2, 0o P 4 [5(k — p) + 5(k + p)]. With sufficient algebra, it is 
possible to show that (3.53) reduces exactly to the L ^ 2 = 0 limit of (3.52). The key 
is in recognizing that 


D 


N 

k 2 


= ( k x q y - k y q x ) 2 (k' 2 + - k 2 ), 

(k + yy^ k 2 + - i/3(k x + q x ) 


(3.56) 

(3.57) 


With our finding that zonostrophic instability encompasses modulational insta¬ 
bility (and the closely related secondary instability), we can envision future avenues 
for research. For understanding how coherent structures grow, a statistical formalism 
like CE2 may provide a clearer window than the fundamental dynamical equations. 
Indeed, a single eigenmode, which is what the calculations from the fundamental 
dynamical equations use, may not be unstable to coherent structures, and instead a 
more complete spectrum may be required. With CE2, one could investigate how ZI 
depends on the background spectrum, using anywhere from a single eigenmode to a 
full incoherent turbulent spectrum. 

In addition, future work could be done to determine how well zonostrophic in¬ 
stability can reproduce modulational/secondary instability when the eigenmodes are 
not Fourier modes, e.g., with nonperiodic boundary conditions. The work presented 
here assumed a Fourier decomposition was appropriate. 

5 There is a seeming factor of 27 t different from the formula in Bakas and Ioannou (2013c) because 
of the choice of Fourier transform convention. 
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3.4 


Beyond Zonostrophic Instability 


3.4.1 Preliminaries: Analogy Between Zonal Flows and 
Rayleigh-Benard Convection Rolls 

The notion of spontaneous symmetry breaking with respect to zonal flows has been 
discussed before (Farrell and Ioannou 2007, Srinivasan and Young 2012). This section 
will expand on that in discussing the mechanics of the symmetry breaking, as well 
as specific consequences it has for the physics of zonal flows (Parker and Krommes 
2013, 2014). 

An important aspect of zonostrophic instability is that it involves a spontaneous 
symmetry breaking. A spontaneous symmetry breaking occurs when a situation’s 
governing physics are invariant under a symmetry transformation but a physical re¬ 
alization is not invariant under the same transformation. A simple example would 
be a ball moving in a symmetric double-well potential, as in Figure 3.5. The equa¬ 
tions of motion of the ball are invariant to reflection about the center line. But with 
friction the ball must eventually end up in one of the wells, a state which breaks the 
symmetry. 

Another well-known example of spontaneous symmetry breaking is the formation 
of convection rolls in Rayleigh-Benard convection (Busse 1978). A box of fluid, taken 
to be infinite in both horizontal directions and finite in vertical extent, is heated from 
below. At weak heating, the heat is transferred to the cooler top surface solely by 
conduction, and the fluid is motionless. But at sufficiently high heating, buoyancy 
forces overcome the inherent dissipation and the conduction state becomes unstable 
to the formation of convection rolls, as shown schematically in Figure 3.6. The con¬ 
vection rolls are spatially periodic but steady in time. This transition to convection 
is analogous to the generation of zonal flows out of homogeneous turbulence. Like the 
conduction state, homogeneous turbulence is (statistically) uniform in space. And as 
a drive parameter such as the strength of the forcing is varied, that uniform state 
becomes unstable to the formation of a periodic structure. Born out of turbulence 
are spatially periodic, steady-in-time zonal flows, which are analogous to the convec¬ 
tion rolls (see Figure 3.7). More than merely descriptive, this analogy will be made 
mathematically precise in the following section. 

3.4.2 Bifurcation Analysis and the Amplitude Equation 

The existence of zonostrophic instability indicates that a homogeneous equilibrium 
without zonal flow is unstable. Perturbations to this equilibrium grow exponentially, 
with wave number dependencies and growth rates that can be calculated. However, 
the ZI calculation alone does not predict how the system saturates. 

To understand the behavior in the regime of nonlinearly interacting eddies and 
ZFs, we turn to a bifurcation analysis. Near the instability threshold, the distance 
from the threshold serves as a small parameter to facilitate analytic progress. It has 
been demonstrated numerically that the bifurcation is supercritical ( arrell and Ioan¬ 
nou 2007), which we confirm with our analytical calculations in Appendix C. Thus, 
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Figure 3.5: Discrete spontaneous symmetry breaking occurs when a ball moving in a 
symmetric double-well potential must, due to friction, end up in one of the wells. 



Figure 3.6: Convection rolls in Rayleigh-Benard convection break the horizontal trans¬ 
lational symmetry. 



Figure 3.7: Zonal flows on a f3 plane break the north-south (statistical) translational 
symmetry. 
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only lowest-order terms in the bifurcation analysis are needed to provide saturation 
of the instability. 

The bifurcation analysis follows a standard procedure, using a multiscale pertur¬ 
bation analysis, expanded around the threshold (Cross and Greenside 2009, Cross 
and Hohenberg 1993). If the threshold occurs at some critical parameter p c , then a 
normalized parameter can be defined as e = (p — p c ) / p c . If we denote u as the state 
vector relative to the homogeneous equilibrium, i.e., u = {W — Wh,U}, then the 
expansion proceeds as 

u = e 1//2 Wi + eu 2 + e 3 / 2 M 3 + ■ ■ ■ . (3.58) 

At first order, one finds 

Ui — A(y, t)r + c.c., (3.59) 

where c.c. denotes complex conjugate [analytically, we work with the quantities 
W(k x ,k y | y) and U(y), which both must be real]. Here, v,\ is proportional to the 
eigenmode r ~ e tQc ^{SW, 5U} that undergoes bifurcation, and A is its amplitude. The 
amplitude is an envelope that slowly varies in space and time. The slow variation 
represents the effect of the infinity of wave numbers nearby q c that also go unstable 
when e > 0. The goal is to determine A, as then U\ will be fully specified. Here, one 
determines a PDE for A as a solvability condition at third order in the perturbation 
expansion. One eventually finds 

c 0 d t A(y,t) = eciA + c 2 <9|A - c 3 |A| 2 A, (3.60) 

where the c; are the order unity, real, positive constants to be calculated. If c 3 were 
negative then one would have a subcritical bifurcation. Equation (3.60) is referred to 
as the amplitude equation, or sometimes as the real Ginzburg-Landau equation. 

It turns out that in order to understand the qualitative behavior of A, one does 
not need to carry out this calculation of the q explicitly ( Toss and Greenside 2009). 
This is because the translation and reflection symmetries (2.36) constrain the lowest- 
order PDE for A to consist generically of the form in (3.60). For example, as a 
result of the translation symmetry, if A is a solution then so must Ae 10 be for any 6. 
This arises because the phase of A determines the location of the solution in space. 
This symmetry requirement demands that the lowest-order nonlinear term is uniquely 
determined to be |A| 2 A. 

The behavior of (3.60) is universal in the sense that, as long as all of the c* > 0, 
the qualitative behavior does not depend of the value of any of the c t . This can be 
seen because all parameters can be transformed to unity by simple rescaling. The 
rescaling is accomplished by letting t = t'T, y = y'L, and A = A'G. One Ends that 
with T = co/eci, L 2 = c 2 /eci, and G 2 = eci/c 3 , that the resulting equation for A' is 
simply 

d t ,A\y\t') = A' + <9|,A' - |A'| 2 A'. (3.61) 

Even if the qualitative behavior is understood, it is still worthwhile to carry out 
the calculation of the coefficients c t . First, computing these and verifying the results 
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Figure 3.8: Comparing showing agreement between numerical solution (blue circles) 
and analytic solution (black line), (a) Compensated growth rate as a function of e at 
q = q c . (b) Growth rate as a function of q at e = 0.01. (c) Compensated zonal flow 
amplitude as a function of e at q — q c . (d) Zonal flow amplitude as a function of q at 
e = 0.0025. In (a) and (c), the compensated growth rate and zonal flow amplitude 
agree with the analytic result at e = 0. The deviation from the lowest order result 
is O(e) for the growth rate and 0(e 1 ^ 2 ) for the zonal flow amplitude. For details, see 
Appendix C. 


numerically provides a concrete check on our overall understanding. Second, the 
perturbation solution may be convenient for certain numerical methods where it is 
useful to start with a good approximation to the true solution. In Appendix C, we 
perform the derivation of (3.60) and obtain expressions for the c*. This computation 
has also been carried out independently (Bakas and Ioannou). To verify our results, 
we compare the analytic growth rate found from (3.60) with that from the exact 
dispersion relation (3.25). Similarly, the analytic ZF amplitude found from (3.60) is 
compared with that from solving the ideal states numerically as in Section 4.1. The 
results are shown in Figure 3.8 and are in excellent agreement. 

With (3.60), the analogy between the zonal flows and the convection rolls in 
Rayleigh-Benard convection is complete. The transition to convection is governed by 
the same class of bifurcation and subject to the amplitude equation. The similarities 
between zonal flows and convection rolls alluded to in the previous section are not 
merely descriptive, but mathematical as well. 

The amplitude equation (3.60) is well understood ( loss and Greenside 2009, 
Cross and Hohenberg 1993, Hoyle 2006), and much of its qualitative behavior is seen 
generically in pattern formation systems. First, with all the parameters c t and e set 
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Figure 3.9: Merging behavior in the amplitude equation (3.60) [ReA(|/, t) is shown]. 


to unity, a steady-state solution exists for any wavenumber within the continuous 
band — 1 < k < 1. To see this, observe that A = ae lk ^ with \a\ 2 = 1 — k 2 is a 
solution. Second, only solutions with k 2 < | are linearly stable; those with k 2 > | 
suffer the Eckhaus instability ( Toss and Greenside 2009). In the Eckhaus instability, 
long-wavelength perturbations grow atop a periodic pattern ( lekhaus 1965, Kramer 
and Zimmermann 1985, Tuckerman and Barkley 1990). This is demonstrated in 
Figure 3.9, where an unstable solution that has been slightly perturbed undergoes 
merging behavior until a stable wave number is reached. Similar merging behavior 
was studied by Janfroi and Young (1999). 

The stability diagram for the amplitude equation is shown in Figure 3.10. The 
neutral curve (N) indicates marginal stability of the A = 0 solution as a function of 
the wavenumber k and control parameter e. The A = 0 solution is unstable to those 
k that are above or inside the neutral curve. At a fixed e > 0, steady-state solutions 
with A ^ 0 exist at any of the k inside the neutral curve. The marginal stability of 
these A ^ 0 solutions is indicated by the Eckhaus curve (E). Inside the E curve is a 
smaller band of wave numbers for which the steady-state solutions are stable. 

Additionally, the amplitude equation is a gradient system, meaning that it can be 
written as 


c 0 d t A(y ) 


5F[A, A*] 
5A*(y) ’ 


(3.62) 


where 

F[A, A*] = J dy ec\AA* + C 2 (dyA)(dyA*) + — c^A 2 A* 2 ^ . (3.63) 


Along solution trajectories F(t) is nonincreasing in time since 




dy 


DA 


dt 


< 0 . 


And F(t) is bounded from below because it can be rewritten in the form 

2 i _2^2 


F\A,A‘] = / iy 


H AA '- e i) + ^ A ' 2 -l e ~i 


1 Ac 1 
~ 2 c 3 


(3.64) 


(3.65) 
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Figure 3.10: Stability diagram for the amplitude equation. The labels ‘stable’ and 
‘unstable’ refer to the nonzero-A steady states. 


where D is the size of the integration domain, which must be selected so that boundary 
terms that arise in integration by parts vanish (eg., D is a periodicity length). Of all 
the solutions A = ae lk A the one with k = 0 gives the smallest value of F. Therefore, 
one might at first think that all initial conditions will tend towards the k — 0 solution. 
However, nonzero k give legitimate steady state solutions, with dA/dt = 0 and hence 
dF/dt = 0. The landscape of F in the function space of all possible A is then such 
that there is a stationary value for each allowed k. Around that stationary value, 
F must be locally flat in one “dimension” corresponding to infinitesimal translation 
and locally increasing in others, but not decreasing since it is a stable equilibrium. 
The k = 0 solution gives a global minimum of F. But even though k — 0 may seem 
to be preferred, this does not guarantee that it is dynamically preferred. Simulations 
with periodic boundary conditions can clearly find nonzero k as the steady state 
solution, as seen in Figure 3.9. This behavior, and more generally the distribution 
of final wavenumbers, has been thoroughly investigated in simulations of the Swift- 
Hohenberg equation, which is also a gradient system (Schober et al. 1986). However, 
in any realistic system, small amounts of noise are present, which perhaps has an 
effect in pushing a physical system towards the minimum of F. It should also be 
noted that even though the amplitude equation is a gradient system, pattern-forming 
systems far from threshold are not in general gradient systems. 

The CE2 system is described by this bifurcation and so near the threshold, and 
more generally, it exhibits solutions existing with a range of zonal flow wave numbers, 
with a certain stability region. In Chapter 4, we numerically calculate the equilibria 
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and stability of nonlinearly interacting turbulence and zonal flows directly from the 
CE2 equations. 
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Chapter 4 

Numerical Calculation of Ideal 
States 


In this chapter, we study the steady-state solutions of the CE2 system (2.21) nu¬ 
merically. As we have learned in Chapter 3, the CE2 system has the mathematical 
structure of pattern formation. This means that there are multiple solutions to the 
equations with differing zonal flow wavelengths and there is an interesting global 
stability behavior of these different wavelengths. Using established techniques from 
the field of pattern formation, in Section 4.1 we compute the nonlinear steady-state 
solutions of CE2 to find self-consistent equilibria of interacting zonal flows and tur¬ 
bulence. We follow that with a calculation of their linear stability in Section 4.2. 
Finally, in Section 4.3 we perform some preliminary exploration into the problem of 
wavenumber selection of zonal jets. 

In the context of an infinite domain with no boundaries, we refer to the steady- 
state solutions as ideal states. Let q denote the fundamental ZF wavenumber of an 
ideal state. For a given q, we solve the time-independent form of (2.21) directly. Our 
approach, which does not involve time evolution, differs from conventional numerical 
studies of turbulence. Time-evolving simulations yield physically relevant, stable 
solutions. In the vast majority of studies these are the solutions one is interested in. 
But when one is interested in the nonlinear dynamics of a system as a whole, one 
often needs to understand the unstable solutions as well. This approach, as well as 
the numerical methods we employ, was successfully used to study convection rolls in 
Rayleigh-Benard convection (Busse 1978). As discussed in Section 3.4, zonal flows 
are mathematically analogous to convection rolls. It is therefore appropriate to use 
the proven techniques on our problem. 

Since we are able to select q and determine the ZF wavelength 2ti/ q directly, 
this method differs from finite-spatial-domain techniques ( arrell and Ioannou 2003, 
2007, Tobias and Marston 2013). Within a finite spatial domain, the wavenumbers 
takes discrete values. The dominant ZF mode is not preselected and is typically 
not the lowest mode because the system evolves self-consistently to find a solution. 
Our infinite-domain technique, which allows the selection of the ZF wavelength, is 
advantageous for understanding the global dynamics. We can solve for both stable 
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and unstable solutions as we continuously vary q and therefore can easily determine 
stability boundaries. 

There is an important side effect to our choice of the dominant ZF wavenumber. 
If the dominant ZF wavenumber does not occupy the lowest mode, subharmonics 
of the dominant mode can be excited. Our method involves setting the lowest ZF 
mode to be the dominant one. This requires fewer resolved modes. It also excludes 
subharmonics (although that is not a limitation in principle). 

4.1 Calculation of Ideal State Equilibrium 

In this section, we compute directly the ideal states of the CE2 system (2.21). Using 
Galerkin projection, we derive a set of nonlinear algebraic equations in a form suitable 
for numerical implementation. The derivation follows procedures first established for 
Rayleigh-Benard convection rolls (Busse 1967, Busse and Clever 1979, Clever and 
Busse 1974, Cross and Greenside 2009, Newell et al. 1990). 

4.1.1 Derivation of Formulas for Numerical Implementation 

The stationary equations for the eddies and zonal flow are 

- (U + - U_)d x w + (u" + - u"_) (V + ^ d x y 

+ [2(3 - (u" + + u”_)]dyd y d x y + F(x, y)-2(fi + uD h )W = 0, (4.1a) 

- [n + v(-l) h d\f]lU{y) - dyd y d x ^( 0,0 [ y) = 0. (4.1b) 

The Galerkin approach begins by expanding the ideal state in some suitable basis 
functions. Appropriate basis functions give rapid convergence as one includes more 
terms, so that one does not have to keep an impractical number of terms. We represent 
an ideal state using a Fourier-Galerkin series with coefficients to be determined, 

MNP 

W(x,y\y)= Y, Z Y, W ^ imax e inby e ipq y 

m=—M n=—Np=—P 
P 

U(y) = V i/y™, 

P=—P 

Here, q is the assumed basic wavenumber of the zonal flow, giving a periodicity 2ir/q. 
While the periodicity in y is desired, the periodicity in x and y is artificial. The 
correlation function W should decay smoothly to zero as x,y — * oo. Therefore, ap¬ 
propriate basis functions in x, y would formally decay at infinity, not be periodic. 
One example of such a basis set would be the Hermite functions. However, we use 
the Fourier basis because of its supreme convenience. Thus, a and b, unlike q, are 
numerical parameters. They represent the spectral resolution of the correlation func- 


(4.2a) 

(4.2b) 
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tion and should be small enough to obtain an accurate solution. Alternatively, the 
box dimensions 2tt/ a and 2n/b must be sufficiently large. 

Because the CE2 equations have translational symmetry in y, there is an infinite 
number of solutions, all equivalent, corresponding to displacements along y. In order 
to obtain a well-posed numerical problem, one must restrict the set of solutions. To 
this end, we again look to the symmetries (2.36). The CE2 symmetries allow us to 
seek a solution for which 


W(x,y | y) = W(-x,-y \ y) = W(x,-y \ -y) = W(-x,y \ -y), (4.3a) 

U{y) = U(—y). (4.3b) 

In other words, we choose the origin of y such that the reflection symmetries hold 
for the solution itself. We find that such solutions do exist. It turns out that this 
restriction does not uniquely specify the solution, as shifting a solution by a half 
wavelength 5y = 7T /q yields a distinct but equivalent solution. Still, this restriction 
is sufficient to make the problem well-posed numerically. Put another way, the above 
condition acts as a way to single out solutions from a family by constraining the 
phase, in lieu of any other boundary conditions. In order for the above symmetries 
to exist in a solution, we also require the forcing to satisfy 


F(x, y) = F(-x, -y) = F(x, -y) = F(-x, y). (4.4) 


Aside from the previous statements, there is also no guarantee that there is a 
unique solution. Indeed, in the zonostrophically unstable regime, once the above 
ansatz with a specific q has been substituted, there are at least two solutions: the 
equilibrium with zonal flows, and the unstable homogeneous solution without zonal 
flows. In some instances we also find other unstable solutions, which may be artifacts 
of the numerical discretization and could be unphysical. 

The constraints in (4.3), along with the conditions that U(y) and W(x,y \ y) are 
real, force U p to be real, and 


= W m,n,-P’ ( 4 ‘ 5a ) 

U- p . (4.5b) 

Furthermore, we take Uq = 0, as that would merely represent a static uniform 
velocity. Equation (4.5a) becomes easier to understand when separated into real and 
imaginary parts (which is done anyway for numerical implementation). Let 


hhmnp — ^- m ,n,p 


U p = 


kbmnp — Fr rlnp 


iF r 


mnp') 


(4.6) 


where E mnp and F mnp are real. Then the symmetries require that 


Emnp E— mn p E rn _ n p E m ^ n _p, 

Emnp E— mn p Fm,—n,p Em,n,—p’ 


(4.7) 

(4.8) 
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Note that F mnp = 0 if any of m, n, or p are zero. The significance of these symmetries 
may be even clearer when the solution is expressed with sines and cosines rather than 
exponentials. One has 


MNP 

W(x,y\y) = £ [ Emnp cos (max) cos (nby) cos (pqy) 

m,n,p =0 

+ Fmn P sin(max) sin (nby) sin (pqy) , 

p 

U(y) = 5 ~2U P cos(pqy ). 
p= i 


(4.9a) 

(4.9b) 


For deriving the nonlinear algebraic equations, the exponential form is much more 
convenient than the sine and cosine form. 

Let us count the number of independent coefficients. For the E mnp , there are 
m = 0,..., M, n = 0 ,N, and p = 0,..., P, for a total of (M+1)(N+1)(P+1) 
coefficients. For the F mnp , there are m = 1,..., M, n = 1,..., N, and p — 1,..., P, 
for a total of MNP coefficients. For the U p , there are P independent coefficients. 
This gives a total of (M + 1)(N + 1)(P + 1) + MNP + P independent, real, coefficients. 

Since ^(x,y \ y) is related to W(x,y \ y), we also write 


V{x,y \y) = C ‘ 


mnp'- 


D imax ginby ^ipqy 


mnp 


and let 


Crnnp — Fr mnp T iH mnp . 


From (2.27), we find the C mnp and W mnp are related by 


bFmnp — ( k + kyky + — k^ | (k — k v ky + — k£ ) (7, 


Vy.Vy 


y J ^mnp 


+ k, T + ( k v + — k 


cL ' n x 


_ 7-2 7-2 

llyll '*-'rnn P ? 


vy ' 2 v 


L d + k x + ( k. 


h 

y 2 y 


a 


mnp 


(4.10) 

(4.11) 

(4.12) 

(4.13) 

(4.14) 


with identical relations between the E mnp and G mnp and between the F mnp and H mnp . 
We have used k x = ma, k y = nb, k y = pq , and defined 


h 2 + = h 2 + + L~ d l 
hl = h 2 _ + Lf, 

h 2 + = kl + i^ky + -k^j , 
h_ = k x + ^ 'k y — — k^j . 


(4.15) 

(4.16) 

(4.17) 

(4.18) 
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We obtain a system of nonlinear algebraic equations for the coefficients W mnp and 
U p by substituting the Galerkin series (4.2) into the steady-state CE2 equations (4.1) 
and projecting onto the basis functions. To demonstrate the projection for (4.1a), let 


i — Jmax inby ipqy 

Ymnp — ° 0 0 

We project (4.1a) onto (p rst by operating with 



(4.19) 


(4.20) 


For instance, the term —(U + — UU)d x W projects to llH p ' mnp U p ’W mnp , where repeated 
indices are summed over, / ‘] p , mnp = — ima8 mp () p ' +p - tfi (<7 + - cr_), cr± = smc(a±n/b), 
and a± = nb — sb ± \p> q. The other terms of (2.21a), as well as (2.21b), are handled 
similarly. In total, we generate as many equations as there are coefficients. 

Appendix D provides the full details of the projection. We summarize the results 
here. It will be convenient to use a shorthand notation where 


k x = ma, (4-21) 

k y = nb, (4.22) 

k 2 = k 2 + Lf, (4.23) 

ky,u = p'q , (4-24) 

ky.w = pq, (4.25) 

ky,u = k^u + azpL d 2 , (4.26) 

h\ = k 2 x + ( k y ± \kyf , (4.27) 

a± = sine , (4.28) 

a± = nb — sb ± pq/ 2. (4.29) 


Using the complex coefficients W mnp , the nonlinear algebraic equations after pro¬ 
jection take the form: 


0 

0 


T (1) TT TIT I JK*) TT C 1 4- T f '° ) c 

1 rstp l mnp u P l rv mnp “T 1 rstp'mnp' J p’^mnp ~r 1 r stmnp^'mnp 


r(2) 


r( 3 ) 


I r( 4 ) TT rr 

' 1 rstp'mnp' J P' ''-’imp 

r(7)r/- r(8) rr 

1 tp’ u P' ' 1 tmnp 11 mnp- 


i d 5 ) _i_ /( 6 ) 
' J-rst ' -*1 


rstmnpWmnpi 


(4.30) 

(4.31) 
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In separate real and imaginary parts, they take the form 


n _ t(1) TT TP i t( 2 ) TT TT _i_ d 3 ^ TT 

u — J rstp'mnp u P ,r mnp "T J rs tp'mnp u P' n mnp ~r orstmnp^mnp 

I t( 4 ) TIM _i_ p( 5 ) _i_ p( 6 ) TP 

' J rstp'mnp u P 111 mnp “r “r °rstmnp^mnpi 


D — A'h 1 ) rr, zt 1 i_ zv( 2 ) rr ,n i zvw n 

O — 1 ^rstp'mnp u P : ^mnp T 1 ^ rs tp'mnp u P'^mnp T r stmnp^* mnp 

I ry( 4 ) 7 r j ,n \ zy( 6 ) F 

“r - fv rstp'mnp u P ,yj mnp ~r i\ rstmnp ± m np , 

a r(7) jT 4- P^ 8 ) Pf 

u J tp' “T 1 1mnp 11 mnp- 


'(3) 


(4.32) 

(4.33) 

(4.34) 


In practice, some of the sums are trivial, and a more convenient form is as follows: 


0 

0 

0 


— t(X) TT F 4- d 2 ^ TT TT _l TT 

~ J rstp'np u p’ r rnp T J rstp'np u P' 11 rnp ~r < J rst 11 rst. 

4- d 4 ) TT TT 4_ d®) 4- p(®) F 

“r J rstp'np u P ,rLrn P ' J rst ' J rst I1 'rsti 

— Zvd) TT , TP -L P<d 2 ) TT .C 4_ zyl 3 )/^ 
t — ly rstp'np u P ,IZ/ rnp ~r r stp'np u P l(J rnp “r l\ rst Ksr rst 

\ zy( 4 ) TT ,C 4- py( 6 ) F 

“ r IY rstp'rnp u P’^rnp T I'-rst 1 rsti 


(4.35) 

(4.36) 

(4.37) 


In these expressions, for J rst and /h r . st there are implicit sums only over p',n,p, but 
no sum over r, s, t. For L p , there are implicit sums over m, n, but not over p. 

In the above expressions, 


d 1 ) _ 

rstp'np rstp 

t( 2 ) _ _ zv( 2 ) — 

u rstp'np v rstp'np 


J rstp'np ^rstp'np k x {&+ & —)$p' +p—t,0i 

kxky : U (°"+ cr -) ( ^p'+p—1,05 

Jrstp'np = —^rstp'np = ky U k x k y k y ^w{ IJ + + O'-^p'+p-^Oj 

where here k x = ra, k y = nb, and 

/( 3 ) — _7y( 3 ) — k k- tt- 
**-rst ^ j b'r yj x ' v y' v y,W i 

42 = K$ = - [2m + ■'(ft? + fc 2 - )]. 

where here k x = ra, k y = s6, k y w = tq. We also have 

I<P = -(p+ *!&)%'„/% 

T (®) — _i. p. u_ 

1 mnp — 


(4.38) 

(4.39) 

(4.40) 

(4.41) 

(4.42) 


(4.43) 

(4.44) 


where here k x = mo, k y = nb, and % = k y jj = pq. 

For the J rst , we have r = 0,..., M, s = 0,..., N, t = 0,..., P. For the K rst , we 
have r — 1,..., M, s = 1,..., IV, t = 1,..., P. For the L p , we have p — 1,,....., P. 
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Schematically, we have the vector of independent coefficients 


x = 



(4.45) 


and the residual vector 

( Jrst \ 

Krst\ 


(4.46) 


We want to solve the system of equations /(x) = 0. 

The system of nonlinear algebraic equations is solved with a Newton’s method 
(Kelley 20CK ). The Jacobian matrix is sparse and is easy to specify analytically, as 
described in the following section. We note that because the ZF equation is linear, 
it is possible to eliminate the ZF degrees of freedom analytically. This is avoided, 
however, because the reduction of only P degrees of freedom is negligible and this 
step incurs the major disadvantage of making the Jacobian no longer sparse. 

A Newton’s method requires a good initial guess. An accurate initial guess near 
the instability threshold is provided by the bifurcation calculation described in Sec¬ 
tion 3.4. To find other solutions we employ simple numerical continuation, where the 
solution at one value of a parameter is used as the initial guess for the solution at the 
next value of the parameter. 


4.1.2 Jacobian Matrix 

It is not too difficult to specify the Jacobian matrix. Take the variation of the residual 
/ by varying the coordinates U p , W mnp in (4.46): 


5 Jrst = 4 3 2SH rst + J®8E, 


L rst 

r(i) 


rst 


n x ; tj 

1 d rstp'np' J p' ur rnp 
I STT 77 > 

“I - J rstp'np u{J P ir rnp 


f( 2 ) 

rstp 

r(2) 


+ J y r2p'npUp'8H rn p 


+ , 5UwH rnv 

1 rstp'np P rnp 


r (4) 

rstp 

K 4 ) 


+ J?slp'npUp'8H rnp 


1 J v±) STJ M 

' J rstp'np UU P ,rL rnp j 


Sl<r„ = K<$6G r x + K%6F r „ 


K l Xn/UrE rnr 


42) 

rstp 

^(2) 


+ K:t P 'n P U P '8G rnp 


+ KZp'np^Up'Gmp 


'(4) 

rstp 

44 ) 


+ K:t P >n P U P '8G rnp 


+ K7tp>nptUp'G rnp , 


5L V = lPsu v + 


(4.47) 


(4.48) 

(4.49) 


Then 


$ Jrst 

5f x (5x) = | 5K rst 
5L r , 


(4.50) 


This gives the Jacobian-vector product at a given point x acting on a vector Sx. 
Implementing this product is virtually identical to implementing the residual vector 
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/ itself. With a little bit of work, one can easily extract the actual Jacobian matrix 
itself, 

Aij = dxj' (4 ' 51) 

The Jacobian matrix is sparse and should be represented as such. When calculating 
the matrix coefficients, one needs to remember to convert the terms 5G —» SE and 
SH ->■ 5F. 

4.1.3 Results 

An example of an equilibrium with p = 0.08 and q = 0.5 is shown in Figure 4.1. 
In the top left is shown the zonal flow velocity U and the strength of turbulent 
fluctuations (measured by the local enstrophy density W ) as a function of y. In the 
top right is a plot of the spectral content of the zonal flow, in both linear and log 
scale. For these parameters, most of the ZF energy resides in the first two harmonics. 
In the middle row is the spectral content of the correlation function W, for the 
homogeneous part W(k x , k y [ p — 0) and the first two harmonics, W(k x , k y \ p — 1) 
and W(k x ,k y | p — 2), of the inhomogeneous part. The external forcing is a thin 
ring in k space around k — 1 that drives only the p = 0 component. The nonlinear 
interactions with zonal flow act to induce a rich structure in the spectral content. In 
the bottom row is W as a function of the real space variables x, y, at several values 
of y. At y — 7t/2 q where the ZF shear is strong, the correlation function in real space 
W(x,y) is distorted compared to its more regular pattern at y = 0 and y = n/q where 
the shear is weak. 

Figure 4.2 shows the ZF amplitude coefficients U p as functions of q at p = 0.21 
and p = 0.19. Near the instability threshold, ideal states exist at all q for which the 
homogeneous equilibrium is zonostrophically unstable [between the two lines labeled 
N in Figure 4.2(a)]. At fixed p, as q approaches the neutral curve boundary (N), the 
zonal flow amplitude falls to zero and the turbulence becomes homogeneous. 

Farther from threshold, there is a region of q where the ideal state solution disap¬ 
pears [between the lines N and D in Figure 4.2(b); see also Figure 4.4], This latter 
bifurcation is not well understood, and may be a result of some other kind of insta¬ 
bility. We believe the feature not to be a numerical artifact. It requires P > 1 to 
exist, but adding more harmonics or refining the resolution do not alter its behavior. 
Moreover, it appears robustly when using multiple variations of Newton’s method as 
well as a distinct Levenberg-Marquardt nonlinear-least-squares algorithm. 

The computational method as described above works very well near the threshold 
p c = 0.237 (7 C = 6.02). However, far from the threshold, for p < 0.12 (7 > 14.2), 
the numerical method breaks down. This appears to be related to the existence 
of multiple solutions at a given parameter value, of which some are unphysical or 
unstable. Far from threshold the Newton’s method seems to inevitably get stuck on 
one of these undesirable solutions. A plot of the spectral content of one of these 
solutions is shown in Figure 4.3. At certain k x values, as k y changes there are strong 
oscillations in the correlation function. This problem does not resolve with higher 
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Figure 4.1: Ideal state equilibrium. See text for details. 


resolution. However, the time-evolving simulations previously mentioned do not have 
these problems because the CE2 equations are statistically realizable and will only 
approach physical, stable solutions. 


4.2 Stability of Ideal States 

With the calculations of the ideal states in hand, we now turn to calculating their sta¬ 
bility. Ideal-state stability, which concerns the inhomogeneous equilibria, is distinct 
from zonostrophic instability, which is a property of the homogeneous equilibrium. 
Both types of instabilities can be described within the CE2 formalism. 
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Figure 4.2: Zonal flow amplitude (7i, f/ 2 as a function of ideal state wave number q 
at (a) n = 0.21 (7 = 7.03) and (b) /i = 0.19 (7 = 7.97). In the unshaded region, ideal 
states are stable. The vertical lines correspond to various instabilities which separate 
the regions (see Figure 4.4). Here, Ld = 00. 
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Figure 4.3: Unphysical solution with strange behavior found by the Newton’s method 
in certain regimes. See text for details. 
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4.2.1 Derivation of Formulas for Numerical Implementation 

Suppose there is an equilibrium {W,U}. We consider perturbations about the equi¬ 
librium: 


W(x,y\y,t)=W(x, y\y) + SW(x,y\ y, t ), (4.52a) 

U(y,t) = U(y) + 5U(y,t). (4.52b) 

The CE2 equations linearized about this equilibrium are 

d t 5W = -( 8U+ - 8U_)d x W - (U+ - U_)dJW 

+ (5Ul - 5U"_) v 2 + ^ d x V + {u" + - u"_) (V + -^dpjd x 8* 

- {5U" + + 5U”_)dyd y d x ^ - (U" + + u''_)d y d y dj^ 

+ 2pdyd v d x 8V - 2{fi + uD h )SW, (4.53a) 

dJSU = -[fi + v{-l) h df]l5U - dyd y d x 5^( 0,0 | y,t). (4.53b) 

With our Fourier-Galerkin solutions in Section 4.1, the underlying equilibrium 
is periodic (in every coordinate x, y, y). Therefore, the differential equation for the 
perturbations is linear with periodic coefficients. If we had imposed periodic boundary 
conditions, then the perturbations would have the same periodicity. But since we are 
assuming an infinite domain, more general behavior is possible. The Bloch Theorem 
states that we can expand the perturbations as a Bloch state ( lever and Busse 1974, 
Cross and Greenside 2009): 

5W(x,y \y,t) = e cr( ' Cl ' q ' >t e lCl ' :x '5WQ(x,y \ y), (4.54a) 

SU{y,t) = e a{Q ' q)t e iQ y v 5U Q (y), (4.54b) 

where the eigenvalue er(Q, q) depends on both Q and q, 5 Wq is the Bloch function and 
has the same periodicity as the ideal state, and Q = (Q x , Q y , Q y ) T and X = (x, y, y) T . 
The Bloch wavevector Q can be chosen to live in the first Brillouin zone: 


— \a < Q x < |a, (4.55) 

-\b < Q y < \b , (4.56) 

~\q <Qy< \q- (4.57) 

We can expand 8Wq and SUq in the same basis functions we used for the ideal state: 

MNP 

SW Q (x,y\y)= ^ SW mnp e imax e inby e ipq y, (4.58a) 

m=—M n=—N p=—P 
P 

tiUQ(y) = 6U P eipqv • ( 4 - 58b ) 

p=—P 
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While the periodicity in the y variable is legitimate, the periodicity in x and 
y are artifacts of the use of a Fourier series. In actuality the correlation function 
should decay as x,y —> oo, not be periodic. If one chooses a and b small enough, 
approximating an infinite domain better and better, one sees that Q x and Q y are 
restricted to lie in a smaller domain near zero. At higher resolution, Q x and Q y 
would presumably get close enough to zero as to not matter. (Nonperiodic basis 
functions such as Hermite functions would not lead to a Bloch wavevector.) A separate 
symmetry argument also suggests taking Q x and Q y to be zero. Due to the correlation 
function exchange symmetry, which we continue to enforce in the perturbations, we 
require SW(x,y \ y) = SW(—x, —y \ y). This requirement forces Q x to be either zero 
or t \a, and Q y to be either zero or | b. To see this, consider a function f(x) of only 
one variable, expressed as 


M 

f(x)=e‘ Ql V y m e imx , (4.59) 

m=—M 

where Q can be chosen to lie within — \ < Q < Suppose / obeys the constraint 
f(x) = f(—x). Then, after reindexing one of the sums with m —> —m, the constraint 
leads to 

M M 

y m e imx = e~ 2iQx y-me imx . (4.60) 

m=—M m=—M 

This equation must be satisfied for all x. It will not be satisfied unless Q = 0 and 
y- m = y m , or Q = \ and y m = y- m - 1 - This fact, when taken with the previous 
argument, strongly suggests taking Q x and Q y to be zero, which is what we do. 
Thus, we have 


5W(x,y | y,t) = e at e lQv 5W Q (x,y \ y ), (4.61a) 

5U(y 1 t) = e" t e iQ y6U Q (y), (4.61b) 

and 5Wq and 5Uq are given as in (4.58). This leads to 

SW{x, y | y,t) = e at 3W mnp e imax e inby e i{Q+pq) y, (4.62a) 

mnp 

5U(y , t) = e at 5U p e i{Q+pq)y . (4.62b) 

v 


The procedure next involves a projection and is similar to that used for the calcu¬ 
lation of ideal states. But several of the symmetry restrictions on the ideal states must 
be relaxed for the perturbations. For instance, the reality condition no longer applies. 
We are looking for a complex Bloch eigenvector. If 6W\ is an example eigenvector, 
real solutions are obtained from 


SW = AWFi + c.c., 


(4.63) 
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where A is some complex amplitude. We also cannot require 5U(y ) = 5U(—y ) or 
5W(x,y | y) = 5W(x,—y \ —y). Furthermore, we should allow for nonzero 5U 0 , 
as there is no reason to discard it in general (except for the case Q — 0, in which 
case 5Uq should be taken to vanish). One constraint that we do retain, as mentioned 
previously, is the exchange symmetry, 5W(x,y \ y) = 5W(—x,—y,\ y ), which is a 
symmetry of all correlation functions. The exchange symmetry requires that 


(fkhrrmp — b\\' — m ,—n,p- 


(4.64) 


Because the eigenvectors themselves are complex, it does not seem beneficial to 
decompose 5W mnp or 5U P into real and imaginary parts, so we leave them as com¬ 
plex coefficients. Let us count the number of independent coefficients. We have 
the 5W mnp , m = —M, ..., M, n = —N, ..., N, p = —P,..., P, with the condition 
that 5W mnp = 5W- m - n:P . Therefore, at each p there is a symmetry much like the 
reality condition of a 2D Fourier transform (but does not involve a complex con¬ 
jugation). In implementation, we choose to keep the following: for n = 0, keep 
m = 0 ..., M, and for n — 1 ,... ,N, keep m = — M, ..., M. This gives a total of 
[M + 1 + N(2M + 1)](2P + 1) complex coefficients from the 5W mnp , and 2P + 1 from 
the SU p . However, one could choose a different implementation such as keeping all of 
the n. The projection of the perturbation equations (4.53) onto the basis functions is 
not affected by the particular implementation of which independent coefficients are 
retained. 

Equation (4.53) is projected onto the basis functions in nearly the same way as in 
the ideal state calculation. The projection results in a linear system at each Q for the 
coefficients SW mnp and 5U P ; this determines an eigenvalue problem for a. Appendix E 
provides the full details of the projection. We summarize the results here. 

In order to relate SW mnp and 5^ mnp , observe from (2.27) that 


fiw mnp 


+ I ky + 2^V’ SW 


Ld 2 + K + ( ky - ~ky,SW 




mnp 


h_| - : swk— ) sw^'^ mn pi 


(4.65) 


_2 ^ 

where k x = ma, k y = nb, k V}S w = Q + PQ, and h ± 5w = L^ 2 + k%+ (k y ± \k yiSW ) ■ 
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It will be convenient to use a shorthand notation where 


kx 

= ma , 

(4.66) 

ky 

= nb, 

(4.67) 

ky,U 

= p'q, 

(4.68) 

ky,W 

= pq, 

(4.69) 

ky } SU 

= Q + p’q, 

(4.70) 

ky,6W 

= Q + pq, 

(4.71) 

k 2 

K y,u 

= h'iu + ®zfLj 2 , 

(4.72) 

T? 

K y,6U 

= ky,su + &ZFL d 2 , 

(4.73) 

k 2 

= k 2 + Lf, 

(4.74) 

h 2 

n '±,SW 

/ 1 \ 2 
= k l + [ky ± 2 k y> sw ) ’ 

(4.75) 

V±u 

. (OL±ulI\ 

= smc l 6 b 

(4.76) 

<y±su 

fa±su^\ 

= smc l 6 b 

(4.77) 

a±u 

= nb — sb ± p'q/ 2, 

(4.78) 

®±6U 

= nb — sb± (Q + p'q)/ 2. 

(4.79) 


After projection, the W equation takes the form 


a5W rst I r stp'mnpfiUp'W mnp + Irstp'mnpUp’SWmnp 


r(2) 


I r( 3 ) rrr .t, 

"T" 1 rstp'mnp uu P , ^mnp 


r(5) 

r(7) 


I rw yrr rrj 

“r" 1 rstp'mnp u u P’ ^ nmp 


+ J (4 ] , 

1 rstp'mnp P mnp 
I 7 ( 6 ) rr 

“r" 1 rstp’mnp u P ,u ^mnp 


+ T KI> AvT/ 

~ 1 rstmnp u ^ mnp 


i A 8 ) AM/ 

' 1 rstmnp u v v mnp • 


( 4 . 80 ) 


As in the calculation of the ideal-state equilibrium, it is more convenient in practice 
to give this formula after performing some of the trivial sums. It becomes 

*8Wr« = I%^Up,W rnp + jlllnpUp'SWrnp 

+ Kstp'np^Up'^ rnp + ^Istp'np^P 1 ^ ™p 
+ Irstp'npfiUp'^rnp + ^Istp'np^P'^rnp 


+ I%6* rat + 4f t SW rst , 


( 4 . 81 ) 


where the implicit sums are only over p',n,p. The zonal flow equation is written 

aSUp = T^pS-Vmnp + T? 0) SU pt ( 4 . 82 ) 

where the implicit sums are only over m, n. 
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In the above expressions, 


(4.83) 

(4.84) 

(4.85) 


Kstp'np — 'i ,k x{ ( 7+,5U °'-,<5£/) ( V+p-t,0; 

Kstp'np — ~^ k x{ a +,U ~ a -,u)^p'+p-t,0: 

Kstp'np = ^ k x k y t SU ^ + ^ k y,W^ ( (T +,5C/ — a -,Su)^p'+p-t,0j 

%%/np = ik £y,u (k 2 + \ k lsw^j {o+,u - 0"-,^) V+p-t,o, (4.86) 

r\ _2 

I rstp'np = ~i' k y,8U k x k y k y,w{®+,5U + cr -,<5t/) ( ^p , +p-t,0; (4.87) 

-—//? \ _2 

I rstp'np = ~'i' k y,u k x k y k y,sw{ cr +,u <?—,u)fip , +p—t, 0) (4.88) 

where k x = ra here, and the other notation is as before. Also, 

I^st = —2i/3k x k y ky } sw, (4.89) 

Irst = — [2^ + U {h, 2 +,5W + k- l ,sw)] J (4.90) 

where k x = ra, k y = s6, and = Q + tq here, and the other notation is as before. 

Finally, 

l£l p = ik x k y k VtB w^ L , ( 4 - 91 ) 

k y,su 

? 10) = - (/* + vkyju), (4-92) 

where ky^su — Q + PQ here, and the other notation is as before. 

If we write the perturbation as a vector 

s * = C W mf) ■ <« 3 > 

then we have an eigenvalue equation for cr, 

aSx = A<5x, (4.94) 


where A is the linear matrix at the equilibrium point x. The sums as written above 
give the matrix-vector product. However, one can also extract the matrix itself with¬ 
out too much difficulty. The matrix is sparse and should be represented as such. 
When calculating the matrix coefficients, one needs to remember to convert the terms 

hT ->■ SW. 

Note that there is a different eigenvalue equation for each Q. For determining 
stability, one must solve the eigenvalue problem for every Q in — \q < Q <\q. The 
equilibrium is unstable if for any Q there are any eigenvalues of A with Re a > 0. To 
calculate this efficiently, Arnoldi iterative algorithms seem to be the best approach. 
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Finally, we point out that these equations contain ZI as a special case, for which the 
equilibrium is the homogeneous one and Q takes on the role of the wave number q. 

It is possible to show two symmetries regarding the eigenvalue, which follow from 
the symmetry of the ideal state equilibrium. They can be verified directly in a 
straightforward, if tedious, way. First, for arbitrary Q, suppose ( 5Wmnp, S Ly } ) is 

an eigenvector with eigenvalue o\ Then the vector [SWmlp, is also an eigen¬ 

vector, with eigenvalue <j*, and 

SWW r = SW^l np , (4.95a) 

SUf = SU^’. (4.95b) 

This guarantees that every complex eigenvalue comes in a conjugate pair. Second, 
suppose at some Q that is an eigenvector with eigenvalue a. Then 

when the Bloch wave number is —Q, the vector {SWmnp , SU^ is an eigenvector 
with eigenvalue cr*, and 

= SW^X-p, (4.96a) 

SU { ~ Q) = ir/'y. (4.96b) 

Thus, for determining stability one actually needs to check only 0 < Q < \q because 
the eigenvalues for negative Q are symmetric. 

4.2.2 Results 

The stability diagram is shown in Figure 4.4. To vary 7, the fundamental dimension¬ 
less parameter defined in (2.13), we change p and hold other parameters fixed (at 
(3 — 1, Ld — 00, v = 10 -3 , h = 4). The stable ideal states exist inside of the marginal 
stability curves marked E, Li, and Ri, which represent different instabilities. The 
Eckhaus instability (E) is a long-wavelength universal instability, present even in the 
amplitude equation (3.60). The Li and Ri curves represent the marginal stability 
boundary for novel short-wavelength instabilities. 

The zonal jets are spontaneously generated by ZI for 7 > y c = 6.02. For 6.02 < 
7 < 14.2, the stability curve is consistent with the dominant ZF wave number observed 
in QL simulations. For 7 > 14.2, we could not calculate the stability diagram with 
this approach due to the aforementioned numerical issues of finding the steady state. 

Part of an unstable eigenvector for the Eckhaus instability is shown in Figure 4.5. 
For this figure, q is just outside of the marginal stability curve, so the equilibrium 
is barely Eckhaus-unstable, and 0 is very small. On the left is the p — 1 spectral 
content of the ideal state equilibrium. On the right is the p = 1 spectral content of 
the eigenvector. The two are proportional. The perturbation is a long-wavelength 
modulation with otherwise the same spectral structure as the equilibrium. The Li 
and R± instabilities have not been analyzed in detail; that could be taken up in future 
work. 
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Wave Number q 

Figure 4.4: Stability diagram for the CE2 equations. Above the neutral curve (N), 
the homogeneous turbulent state is zonostrophically unstable. Ideal states are stable 
within the marginal stability curves (circles) E, L l5 and Ri. The circle-points were 
computed by solving for ideal states using the method of Section 4.1, then finding 
the marginal stability boundary using the method of Section 4.2. Also shown is the 
dominant ZF wave number from independent QL simulations (crosses). The crosses 
were computed by performing direct numerical simulation of the QL equations (2.17) 
and determining the dominant ZF wave number. The stability region calculated from 
CE2 is consistent with the ZFs realized in the QL simulation. The stationary ideal 
states vanish to the left of curve D. Discussed in Section 4.3 are the interior curves: 
Rhines wave number (black dashed line), wave number of maximum growth rate for 
zonostrophic instability (blue dotted line), wave numbers of minimum eddy energy 
(red line), minimum total energy (black line), minimum eddy enstrophy (red dotted 
line), and minimum total enstrophy (black dotted line). Here, a = 0.06, b = 0.08, 
M = 20, N = 27, P = 4, and other parameters are given in the text. 7 is varied by 
changing /i while holding the other parameters fixed. 
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Figure 4.5: Left: the W(p = 1) component of the ideal state. Right: the SW(p = 1) 
component of the instability eigenvector. The eigenvector of the Eckhaus instability 
is a long-wavelength modulation with otherwise the same spectral structure as the 
equilibrium. 


4.3 Wavenumber Selection 

As evident from Figure 4.4, we are presented with the theoretical quandary of having 
a wide range of allowed, stable solutions and yet a narrow preferred region where QL 
realizations tend to appear. This is common to pattern-forming systems, and the 
problem of wavenumber selection is difficult ( Iross and Greenside 2009). The Rhines 
wavenumber (1.26) can be estimated by using \U 2 = E and E = e/2/i to give 

k R « /^/y/V 1 / 4 . (4.97) 

This estimate works well in giving the preferred ZF wavenumber. In this section we 
explore what features of the equilibrium might correlate with the preferred wavenum¬ 
ber, in an attempt to achieve a greater understanding of what determines wavenumber 
selection of ZFs (Parker and Krommes 201 ). 

One might naturally inquire as to whether the preferred mode is the fastest grow¬ 
ing mode in the ZI about the homogeneous equilibrium. This does not appear to be 
the case away from the threshold at larger 7 (Farrell and Ioannou 2007, Srinivasan 
and Young 2012), as seen in Figure 4.4. There is, however, a plausible scenario that 
emerges which may explain the merging of jets often observed in the beginning stages 
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of simulations, especially those which initialize everything at low amplitudes. At 
large 7 , it appears that the fastest growing mode may be to the right of the stabil¬ 
ity region. In a simulation, the turbulence quickly comes to a quasi-equilibrium on 
a short time scale and begins to drive the ZF. The growing ZF mode cannot stably 
saturate, for its wavelength is too small to coexist with the turbulence. As the system 
evolves through the subsequent instability to drive the jets toward larger wavelength, 
a space-time visualization such as that in Figure 2.1 displays merging jets. 

Another possibility is that some kind of variational principle applies. The ampli¬ 
tude equation, by which CE2 is governed near threshold, is a gradient system. The 
ideal states of varying wave number q have varying values of the effective free energy. 
However, the minimum of the effective free energy is not necessarily dynamically pre¬ 
ferred (Schober et al. 1986). In any case, away from threshold CE2 is not a gradient 
system and there is no rigorous theoretical basis for expecting variational behavior 
to occur. From a different perspective, variational principles for certain 2D turbulent 
systems have long been discussed theoretically. Some of these principles are based 
on the nonlinearly conserved quadratic quantities, the energy and the enstrophy. For 
instance, in freely decaying turbulence where viscosity provides the dissipation, the 
enstrophy is expected to decay more quickly than the energy. One might expect the 
decaying turbulence to reach a state of minimum enstrophy subject to the constraint 
of constant energy. Other principles exist based on minimum dissipation or maximum 
entropy or entropy production (Majda and Wang 2006). Although these principles 
do not directly apply to the damped, driven CE2 system, they at least motivate a 
numerical exploration to try to discover any correlation between the preferred wave 
numbers and other properties. 

As a simple starting point for our exploration, we examine the energy and enstro¬ 
phy of the ideal states. Plots of the energy and enstrophy, for both the total and 
just the eddies, are shown in Figure 4.6 for /i = 0.15. For each quantity a distinct 
minimum is present. We find at each /i the minimum of all four quantities; the result¬ 
ing curves are shown in Figure 4.4. While the minima of the total energy and total 
enstrophy are consistent with the QL realizations, there is no clear indication that 
either is especially preferred. On the other hand, the accessible regime investigated 
here is not too far from threshold, so this is not in the asymptotic regime of large 7. 

There is no definitive conclusion to draw from these explorations. Determining 
and understanding the length scale of zonal flow remains an important and unsolved 
problem in plasma physics. Future investigations along the line discussed here may 
prove useful. 
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Figure 4.6: The energy and enstrophy, both total and that of just the eddies, for ideal 
states of varying wave number q at /r = 0.15. 
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Chapter 5 

A Pedagogical Closure for 
Homogeneous Statistics 


In this chapter we propose a closure for homogeneous statistics. 1 The material here 
is separate and not immediately connected to the work on zonal flows in the rest 
of this thesis. With some extensions described in the following paragraphs and in 
Chapter 6, it could be connected to zonation and inhomogeneity. 2 But even without 
those extensions, it is interesting in its own right and has pedagogical value for its 
simplicity. So while the work here lays the groundwork for further investigation into 
zonal flows, as it stands, it is simply an interesting venture into a turbulence closure. 

Systematic closures for homogeneous turbulence possess important properties, in¬ 
cluding statistical equilibrium, realizability, and an H-theorem ( larnevale et al. 1981, 
Krommes 2002). Here, we introduce a simple closure which exhibits some of the im¬ 
portant properties in a simple, transparent way. 

Additionally, we discuss another property that has not previously received much 
attention, which is the stability of the closure’s steady-state solutions. From our 
discussion of zonostrophic instability, it is clear that stability of the homogeneous 
equilibrium plays a critical role. In CE2 with external forcing and linear damping, 
this stability is trivial. In the closure we introduce here, it is decidedly nontrivial, 
yet we are still able to prove general statements about stability. Finally, the closure 
is simple enough that its equilibria can be completely characterized. We find that 
there is a unique physical solution, and the necessary and sufficient conditions for its 
existence can be explicitly stated. When the wavenumbers are discretized, there is a 
large number of nonphysical solutions. In some special cases, analytic solutions can 
be found. 

This “toy” closure is not intended to be an accurate portrayal of turbulence. It 
is a stepping stone, like the works of Kraichnan and Spiegel or of Leith in the early 
days of analytic turbulence theory (Kraichnan and Spiegel 1962, Leith 1967), in an 
attempt to understand some piece of the puzzle. Those studies used simple closures 

lr This work is unpublished. 

2 This material might be more logically placed as part of Chapter 6, but that would bog that 
chapter down. It could be put in an appendix, but it was desired that this material not be doomed 
to languish in obscurity. 
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of homogeneous turbulence to try to understand the inertial range in 3D. Here, our 
ultimate goal is to get at the fundamental mathematical structure that underlies the 
bifurcation at which zonal flows are born. 


5.1 Model Closure 


Consider a 1-field turbulent model, 


ipk = +Afk[i>k\, 


(5.1) 


where A k represents linear terms, including drive and dissipation, and J\f k incorporates 
all nonlinear terms. We will consider the two-point correlation function, C k = (V’fcV’fc)- 
As an equation that might represent the dynamics of Ck, our closure is 


C k 


jkCk — \x k CkC + f k C , 


(5.2) 


where y*, = Re A*,, and C is independent of k , has the same dimensions as Ck, and 
is in some way a measure of the turbulent intensity. This form is intended to mimic 
the structure of turbulent damping plus nonlinear forcing, which appears in more 
sophisticated closures like the Direct-Interaction Approximation (DIA) or the Eddy- 
Damped Quasi-Normal Markovian Closure (EDQNM). Thus, we require the fk to be 
positive, and we also find that we will have nice properties if all of the fik are positive 
as well. The fk and /p, are left otherwise totally unspecified. 

If energy conservation among the nonlinear terms is desired, there is a unique 
choice for C which provides this. Suppose the energy in a given mode is 

E k = a k C k , (5.3) 

where a k is a positive-definite weighting. Then the total energy is E — a k C k , and 
the rate of change of energy due to nonlinear terms is 


E NL = 


j>*a| 


NL 


= J>* (-HkC k C + fkC 

k 

= C ( C ^ a k fk - ^2 a kHkC k j . 
V k k / 


Therefore, for E NL = 0, we find 


^k^kCk 

Y2k a kfk 


(5.4) 


(5.5) 
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Upon defining 


(5.6) 


F = ^ ffk/t. 

k 

we write 

C = ^ cr k HkC k . (5.7) 

k 

Thus, C is a weighted sum of the energy in each mode (weighted by //*), with a 
certain normalizing factor. This closure yields only one quadratic quantity conserved 
by nonlinear interactions. More sophisticated closures like the Direct-Interaction 
Approximation can conserve multiple such quadratic quantities ( 'Craichnan 1959). 

5.2 Properties 

It turns out that this simple closure exhibits many desirable properties of statistical 
closures. These properties are: 

1. One can prove certain stability properties of the steady-state solutions. In this 
case, every nonzero equilibrium is linearly stable. 

2. Equipartition solution is possible for statistical equilibrium with no linear terms. 

3. The nonlinear closure terms exhibit an H-Theorem for monotonic increase of 
entropy towards the statistical equilibrium. 

4. The system is statistically realizable. That is, if the system initializes with all 
nonnegative Ck, they stay nonnegative. 

5. The nonlinear equilibria of the system can be completely characterized. It 
can be shown that there is one, and only one, physically allowable equilibrium 
(with all the Ck nonnegative), and a necessary and sufficient condition on the 
parameters for its existence can be derived. 

There are of course deficiencies owing to the simplicity of the treatment of the non¬ 
linear terms. For instance, there is no real mode coupling among triads. Only one 
quadratic quantity can be conserved. The closure also does not account for the effect 
of linear waves. 

5.2.1 Linear Stability 

Clearly, Ck = 0 for all k is one equilibrium. However, if any of the 7 k are positive, 
that equilibrium is unstable. We now show that any and all nonzero equilibria are 
linearly stable. We do this by providing a positive definite functional, quadratic in 
the perturbation, which decays in time according to the linearized system. That is, 
small deviations from the equilibrium must eventually die away. This happens to be 
true even if the equilibrium is nonphysical (meaning the proof goes through even if 
some of the Ck are negative). 
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To prove linear stability of any nonzero equilibrium, first note the “steady-state 
condition” of an equilibrium, 

Ik = 1-ikC - A C 2 . (5.8) 

L-fc 

This relation will be used to eliminate 7 *, later on. 

Linearization 

Linearize (5.2) about an equilibrium C k , by letting C k (t) = C k + 5Cj.it ): 

5Ck — 7 k5Ck — fi k 5C k C — n k C k 5C + 2f k C5C, (5-9) 

where 5C k = d5C k /dt, and 5C = J2 k a kl^k5Ck/F. Now substitute for 7 *, using the 
steady-state condition in (5.8) to obtain 

SC k = -fi k C k 5C - A C 2 5C k + 2 f k C5C. (5.10) 

It will be convenient to write this in terms of w k = 5C k (t)/C k : 3 

w k = -fi k 5C - A c 2 w k + ^ C5C , (5.11) 

where now 5C = Y, k °kHkC k w k /F. 

Quadratic Functional 

Consider the quadratic functional 

Wjt) = h k C k w k (t ) 2 , (5.12) 

k 

where h k is a positive quantity. Observe that W is positive definite. We take 

h k = y (5.13) 

3 Technically, for this to be allowed, none of the equilibrium C k can be exactly zero 
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as the weighting factor. Then, the evolution of W is given by 


W 


Y °^C k w k w k 

k 

- Y y (f*kC k w k 6C + f k C 2 wl - 2 f k Cw k 5c) 

k 



kk k ^kC k W k 5C kT k f k —^2 2 

-6 W u 


+ 


2 CFC°— 


Y (Tk -f' kWk 


c 2 


Y (T kfkW 



SC - ^ Y °kfkW k 


2 —2 
_ 

F 2 


^ ^ kkj 17 1; .fj f k k ' j ^ ' k 


jk 


c 2 


^ ^ k^k fk k ' 


SC - ^ Y ak fk Wk ] ~T^[ F Y ak -f' kWk ~ Y a 3 <T kfjfkW j wA . (5.14) 


J 


J k 


In the last equality, the terms in the second set of parentheses can be written 


F Y CT kfkU’l - Y VjVkfjfkWjWk = Y( a j a kfjfkwl - (TjCkkfjfkWjWk) 

k jk jk 

= 2 Y Wkfjfk{2wl - 2 WjW k ) 

jk 

= \Y 0 j°kfjfk{w 2 k - 2WjW k + Wj) 
jk 

= \Y a j a kfjfk{wj - u>k) 2 , (5.15) 

jk 

where to get to the third line, we have swapped indices j k in one of the terms 
(with symmetric combinations in front not changing). 

Thus, we find 


IT 



C 

F 


^ ^ kkkfk^Vk 


c 2 

2 F 2 


^ ^ kkj kk k f j f k ( Wj 
jk 


u>k) 2 , 


(5.16) 


and hence, IT < 0. And IT = 0 only when the perturbation w k = 0. 

Note that normalizing the perturbation 5C k {t) by the equilibrium value C k was not 
merely convenient. It also served to eliminate consideration of the zero equilibrium, 
which is not stable if any of the y*, are positive. What we have shown is that all 
nonzero equilibria are stable. 
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External Forcing 

External forcing could be added to the system, where instead of (5.2), one might have 

Cfc = 7 kCk — kkCkC + f k C + J-fc, (5-17) 

where T k is some known external quantity and is nonnegative. In this case, C k = 0 
for all k is obviously no longer a solution. However, in a minor modification of the 
above proof, one still reaches the conclusion that any nonzero equilibrium is linearly 
stable. 

5.2.2 Statistical Equilibrium 

If we consider (5.2) with all the 7 *. = 0 and just the nonlinear closure terms, we find 
that a steady state is given by 


0 = -fi k C k C + f k c\ (5.18) 

or 

C k = —A, (5.19) 

kk 

where A is some constant and is equal to C(t 00 ). Note that if f k and kk are 
proportional, that is, if f k /kk is independent of k, then we in fact have an equipar- 
tition statistical equilibrium. Also, if any of the [i k are negative, then the statistical 
equilibrium predicts a negative C k and is unphysical. 

By using the fact that energy is conserved, we can determine A from the initial 
conditions. We have 


E(t = 0) = a k C k (t = 0), 

k 

E(t —>■ 00 ) = CT k C k (t — * 00 ) = A °kJ k 
k k V k 

and thus, by conservation of energy, 

A = E k a kC k (t = 0) 

Efc ®kfk/kk 


(5.20) 

(5.21) 


(5.22) 


5.2.3 Entropy and H-Theorem 

Carnevale et al. (1981) have shown that in many second order closures, it is useful 
to think of the entropy as S = '}2 k In C k . Here, it turns out we need to use a slightly 
different definition, namely, the entropy-like quantity 

s = — hl ° k - ( 5 - 23 ) 

k V k 


81 





(We do not prove that this is equivalent to an entropy for this model, but it has 
similar behavior.) This form, like that of Carnevale et al., is scale-independent. That 
is, if C k —> h k C k , the h k terms affect only the absolute entropy but not changes in 
entropy because the (cr k f k /n k ) In h k factors are simply constants. 

The evolution of the entropy is given by 


* = £ 
k 

= £ 


O'kfk 1 
Fk C k 

®kfklk 

Fk 


c k 


£ 


Okfk 

Fk 


~FkC + 



(5.24) 


The contribution from the nonlinear piece can be written 
O'kfk 


£ 


Fk 


^c+^-c 2 ) = c^2(-a k f k + 


vkfk ^ 


C 


And hence, 


C 

F 

C 


Fk C k 


£ OjOkSjSk + £ Ff=r a H i i C i 


£ 

jk 


jk 


O jO k 


jk 


FkC k 


Fine., 

Fk C k 


fjfk 


C y". OjOkfj f k 

2 F FjFkCjCk 

JK 

C OjOkfj fk 
2 F , k FjFkCjC k 

C \ ^ <J j a kfjfk 
2F jk FjFkCjC k 


f! 


o FjFk r< n 
Ffk j k 

J J J r' j 


FjC] _ 2 HFk c , Ck + F-lCj '] 


r 

J 3 

FjCj 

fj 


fjfk 
FkC k \ 2 

~17) ' 


fk J 

(5.25) 


S = 


ST' Qkjkjh C \ ^ OjOkfjfk ( FjCj Fk.Ck \~ (^ 26 ) 

V ^ 2F * FjFkCjCk v fj fk J ' J 


First, notice that there is an H-Theorem. If all the 7 k are zero, then S is positive 
dehnite, and increases monotonically towards the statistical equilibrium state where 
Ck OC fk/Fk■ 

With the 7 k included, then since the contribution from the nonlinear piece is 
always positive, a necessary condition for a steady state to be reached is that the 
contribution from the linear piece be negative: 

OkfkFk q 
Fk 



(5.27) 
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We will see later from a detailed analysis of the solutions that this is also a sufficient 
condition for an equilibrium to exist. 


5.2.4 Realizability 

Proof: Algebraic Assume that the ODE is solved as an initial-value problem, and 

all the Ck are initially nonnegative (and hence C is also nonnegative). Then, as the 

system evolves in time according to (5.2), the C k remain nonnegative. For if Ck were 
_2 

ever 0, then Ck = fkC > 0. Thus, the system is realizable. 

This argument does not easily generalize to models with more than one held, 
however, because if C k = 0 , then linear coupling to a different held C k may still in 
principle cause C k to possible become negative. 


Proof: Langevin Equation Alternatively, one can prove realizability by provid¬ 
ing a Langevin equation which has the same statistics as (5.2). The derivation is 
somewhat similar to that used for the standard Markovian closures in that the terms 
in the Langevin equation depend on the statistics of the solution (but there is no 
analog to the triad interaction time here). Consider the random equation 


where 


^ = 2L 


7fc'0fc ~ rjk{t)il>k 


+ o,k, 


Vk(t) = h kC(t ) = y WAaA, 

k 

a k {t) = w{t)C{t)y/J k . 


(5.28) 


(5.29) 

(5.30) 


We assume the Ck (and C) in the random equation are “known”, or at least nonran¬ 
dom and independent of -0*,. This comes from assuming that in the limit of a large 
ensemble, each individual - 0 k contributes only infinitesimally to the statistics to give 
Ck- Also, w(t) is Gaussian white noise with zero mean and 

(■ w(t)w(s )) = 5(t — s ). (5.31) 

Then, let us compute the second-order statistics for 

C k (t) = (5.32) 

We have 

C k (t) = 2R e(^ k (t)ip* k (t)^ 

= Re(( 7 k - r}k)'$k'&*k) + 2R e(a k (t)^* k (t)) 

= (Re 7 k )C k - (R erj k )C k + 2 R e(a k (t)ijj* k (t)). 
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We assume 7 and r/k to be real. To compute we write down the Green’s 

function solution of (5.28), 


rt 

Mt)=M0)e J o dT ^+ / dse^ dT ^ (T) a k (s), 

Jo 

(5.33) 

where 7 fc = ( 7 k — rjk)/2. We assume ^(0) = 0. Thus, 


««= f isel>^al(s), 

Jo 

(5.34) 

(ak(t)$Z(t)) = 1 dse^ dT ^ k{T \d k (t)d* k {s)), 

Jo 

(5.35) 

(d k (t)d* k (s)) = 5(t- s)f k C 2 . 

(5.36) 


(a k {tWk(t)) = fkC 2 f ds5(t - s)e^ dT ^ 

Jo 

= i hc\ 

Finally, we have 

Ck = ikCk — HkCkC + fkC . (5.37) 

which matches (5.2). 


5.2.5 Nonlinear Equilibria 

Now let’s actually try and solve for the equilibria of the system in (5.2). In equilib¬ 
rium, we have 

0 = (7 k-HkC)C k + f k C 2 . (5.38) 

_ 2 _ 

Since f k C > 0, any equilibrium must have ( 7 *, — [i k C)C k < 0. Therefore for any 

physical equilibrium with all the C k > 0, we must have 

»kC - 7 fc > 0 (5.39) 


for every k. 

From (5.38), we may write 


C k (C) 


fkC 2 
HkC 7 k 


(5.40) 
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which means that if C is known, then Ck for each k is known. The method of solution 
is now to solve for C. Multiply (5.40) by 07 / 2 *, to obtain 


Then sum over k to give 


r , °kHkfk 7^2 
^kf^k^k -p- w • 

Hk^ ~ 7 k 


^ ^ ^kH'kfk 


(5.41) 


k'kC — 7 fc 


-=2 


Dividing through by C (assuming we don’t want the trivial solution with all the 
Ck = 0 ) gives 

F E W.fo ■ (5.42) 


c ^ C — 7 */ Hk 

Substituting in the form of F and combining terms gives 

@kfk'~ik/k'k 
k C 7 k/k'k 


° = E 


(5.43) 


This equation completely describes the nonlinear solutions. One finds the solution 
C, and then computes Ck from (5.40). The above equation has more than one possible 
solution. If there are N modes in the system, then k — 1,..., N. If one multiplies 
through by all the denominators to obtain a polynomial equation in C, one finds a 
degree N — 1 polynomial, and hence N — 1 possible solutions to this equation (adding 
back in the trivial solution C = 0 gives a total of N). To determine how many roots 
are real or complex, we must do something different. Note that statistical equilibrium 
is retained here by noting that for 7 k = 0 , any C is an allowable solution. 

A graphical approach is fruitful. Let x = C, and 27 = 7 &//!&. The problem is 
then restated as solving for the values of x which satisfy 


E 


&kfk%k 
X ~ X k 


0 . 


(5.44) 


Graphically, this amounts to finding where the function on the LHS crosses the x 
axis. For simplicity, assume that all of the 27 are distinct and nonzero. Some of 
the Xk are positive and some are negative, corresponding to positive and negative 
7 *;. Then the graph of this function has N vertical asymptotes, one at each 27 . For 
instance, with N = 4, the graph may look something like that depicted in Figure 5.1. 
Supposing that m of the 27 are negative and n are positive, then because of the 
vertical asymptotes there are guaranteed to be at least (m — 1) + (n — 1) = N — 2 real 
roots. Since we have proven that there are only N — 1 possible solutions, then only 
one other solution can exist, and so it must be real. It is not immediately obvious 
whether the remaining root will occur to the left of all the 27 , to the right of all the 
Xk, or between the two 27 asymptotes surrounding x — 0. We shall discover that 
there is never a root between the two 27 asymptotes surrounding x = 0 , and that the 
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Figure 5.1: Schematic graph of (5.44) with N = 4, with zeros highlighted. For 
&kfk%k ^ o. 


remaining root always occurs to the right of all the Xk, or to the left of all the x fc , 
depending on a certain criterion. 

To determine if a root exists to the right of all the Xk, we simply want to know 
if the function crosses zero. Since the function is monotonically decaying, we can 
answer that question by determining whether it is positive or negative in the x —> oo 
limit. Taking the large |x| limit of (5.44), we find the LHS goes as 


1 

x 


N 

^ ^ ^h•t'k ■ 
k =1 


(5.45) 


If this is negative at large x, there must be a root to the right of all the Xk- This 
condition is given by 

N N „ 

Y (?kfkx k = Yj ° k fc7fc < ( 5 - 46 ) 

k =1 k =1 ^ k 

On the other hand, if that quantity is negative at large negative x, there must be a 
root to the left of all the Xk- This condition is given by the opposite, 

N N f 

= E —> 0- 

tt si w 

Thus, assuming that sum is not exactly equal to zero, one of these two conditions 
must be true, and the final root occurs either to the left or to the right of all the Xk, 
and not between the two Xk asymptotes surrounding x — 0. 
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Now that we have a picture of where the solutions for x or C are, let’s determine 
which solutions are physically allowable. From (5.40), we see that for Ck to be 
positive, the quantity HkC — 7 k must be positive for all k , or equivalently, 

x — Xk > 0, for all k. (5-47) 

If x > Xk for all Xk, then graphically this corresponds to the root of (5.44) being to 
the right of all the vertical asymptotes. We have already shown that the condition 
in (5.46) is the necessary and sufficient condition for such a root to exist. If it does 
exist, it is the unique physical solution. The other solutions correspond to some of 
the Ck being negative. The number of negative Ck correspond to how many of the Xk 
are to the right of the root x. Incidentally, even the unphysical equilibria are linearly 
stable. If the condition in (5.46) is not satisfied, then the system does not saturate; 
it blows up. One may think of the 7 *, as being too large in this case. 

It is remarkable that the equilibria of (5.2) can be fully characterized. With 
a statistical closure, one is primarily interested in the actual steady state, not the 
transient evolution to the steady state. Numerically evolving the statistical closure 
model, like (5.2), can be time-consuming, mainly because a small timestep is required 
to ensure both accuracy and stability. It might be advantageous to write down the 
equation for steady state and find a way to directly compute the solutions. However, 
in general certain difficulties arise when attempting this route. In particular, 

• A nonlinear solution may be unphysical (see, e.g., Section 4.1). Obviously, one 
can discard any unphysical solution that is found, but there is still the problem 
of how to find a physical solution. The time evolution method does not have 
this problem, so long as the model is realizable. 

• A nonlinear solution may be (linearly) unstable. One has no way of knowing 
(without further computation) whether the solution that was found is stable to 
small perturbations, and an unstable equilibrium has no relevance. A separate 
stability calculation can be done (see, e.g., Section 4.2), but this may be difficult. 
The time evolution method does not have this problem, since any equilibrium 
found that way must be stable. 

• Not knowing how many physical solutions may exist. The time evolution 
method also has this problem—an equilibrium may be found, but it is not 
known if others exist. 

The closure used here is simple enough that it can be analyzed in sufficient detail to 
overcome all three of these difficulties. 

I 11 practice, it is not advisable to try to numerically solve the polynomial form of 
the equation for C. This is because finding the roots of high degree polynomials is 
an ill-conditioned problem (Wilkinson 1994). Since roundoff error is inevitable, large 
errors can result, including finding complex roots even though they should all be real. 
Instead, using a standard nonlinear root finder on the form in (5.44) is preferable, 
especially because it is known that there is one and only one zero in the domain 
(max(xfc), 00 ). 
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5.2.6 Continuum 


The same results hold when using a continuum rather than discrete modes. The 
results are summarized below: 




C{k) = 7 (k)C(k) - fi(k)C(k)C + f(k)C , 
E(k) = a{k)C(k), 

E — J dkE(k), 

C — — [ dka(k)/j,(k)C(k ) , 


F= dka(k)f(k). 


(5.48) 

(5.49) 

(5.50) 

(5.51) 

(5.52) 


Then the nonlinear terms conserve the total energy E. Here, k can be a vector. 
Equilibria are obtained from 

o = [ 7 (jfe) - fi(k)C\ C(k ) + f(k)C 2 . (5.53) 


Any physical equilibrium must have fi{k)C — 7 (k) > 0 for all k. Divide through by 
n(k)C — 7 (k) to obtain 


C(k) 


n_k)c 2 

n{k)C — j(k) 


(5.54) 


The above step is only valid for C ^ max (7 (k)/fi(k)) for any k. This is satisfied if 
C > max(y (k)/n(k)). Multiply by a(k)fj,(k) and integrate over k to obtain 


FC 


This can be rearranged as 


F 

U 




a{k)n{k)f{k) 
fi>(k)C — j(k) 


°(k)f(k) 

C--f(k)/ii(k)’ 


(5.55) 


(5.56) 


()= f dk °(k)f(k)i(k)/y(k) 

J C — 'y(k)/n(k) 


(5.57) 


In the discrete case, we could use the properties of polynomial equations to prove 
that there was only one possible solution for C > max( 7 fc/ fik)- Here in the continuum 
case, we have not found a proof that a solution is unique, though that does not mean 
a proof does not exist. 

In the discrete case, we also found a criterion that was both necessary and suffi¬ 
cient for a unique solution. In the continuous case, that proof gives only a sufficient 



condition for a solution. To see this, rewrite (5.56) as 


x J x — 7 


(5.58) 


where C = x and 7 = 7 /ji. If the RHS is smaller than the LHS as x —> 00, then 
there is guaranteed to be at least one solution. If we expand the RHS for small x, we 
obtain 

[ dk — Ac = — [ dkaf H—- [ dkafj + --- 
J x - 7 X J X 2 J 

Fir 

= - + - / dk afj + ■ ■ ■ . 

X X 2 J 

Therefore, the RHS is smaller than the LHS above as x —> 00 if 

[ dka(k)f(k) 7(fc) = f dk vWf&yvW . < 0. (5.59) 

Since we have not yet proven the solution is unique, this approach does not show 
that the criterion is a necessary condition. However, entropy considerations in the 
continuous analog to calculations in Section 5.2.3 do prove it is a necessary condition. 
To find an equilibrium, the nonlinear equations to solve are 

or 

f dk dMtMM = o. 

7 x ~ 7 \k) 


(5.60) 

(5.61) 


A Solvable Example 

Here we provide an example that can be integrated directly and solved for x. Suppose 


a(k) = 1 , 

(5.62) 

W = i + V 

(5.63) 

1 -(3k 2 

tW = t "i + /» 2 ’ 

(5.64) 

"Wm+'V 

(5.65) 

7 (k) = —(1 -/3/c 2 ), 

Ro 

(5.66) 


where fo,^o,7o,/d > 0 and let k extend from 0 to 00 (or from —00 to 00 ; it doesn’t 
change the result). This y(fc) peaks at k = 0 and is one-to-one. The max of 7 is 










7o/ho, so x > 7o//io is required. The nonlinear equation for a steady state is 



fo (7o/ho) (1 ~/3k 2 ) 

1 + (3k 2 x - ( 70 /ho) (1 - /3k 2 ) 


0 . 


Let 2 = x/io /70 — 1 > 0. Then 



1 - fik 2 1 
1 + f3k 2 z + /3k 2 


0 . 


Change integration variables to y = /3k 2 , to obtain 



1 1 - y 

yl/2 +!)(?/ + z ) 


0 . 


We use the integrals 


dy 


7T 


1 0 y 1/2 (y + a){y+ b) y/ab(y/a + y/b)’ 

M 2 

dy 


y 


7 r 


{y + a){y + b) y^+v^’ 


valid for a, b > 0. Equation (5.69) becomes 

7r 7r 


y/z( 1 + \fz) 1 + \fz 


(5.67) 

(5.68) 

(5.69) 

(5.70) 

(5.71) 

(5.72) 


or 


1 = Vz, 



l=z, 


■> ^ 

O 

O 

1 

so that 

x = C = 2^. 

ho 

Substituting into C(k), we have 



C(k) = 4 


7o/o 1 
h 0 1 + /3k 2 ' 


(5.73) 

(5.74) 

(5.75) 

(5.76) 


(5.77) 


A slight modification to the previous example with q(/c) = 70(1 — (3k 2 ) is also 
solvable. This 7 (k) gets continuously more negative at large k, like viscosity, instead 
of saturating at a constant negative value. 
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Another Solvable Example—2D Isotropic 

It is also possible to construct an integrable example that is 2D and isotropic. Let 


a(k) = 1, 

(5.78) 


(5.79) 

in l ~P kl 

7W = T0 1 + /»*’ 

(5.80) 

**> = 1 +v 

(5.81) 

l{k) = — (1 - fik 4 ). 

AD 

(5.82) 


Note that if f(k) went like k 2 at large k, the integral for F would not converge. The 
nonlinear equation to solve becomes 



k 1-13k 4 
l + /3k 4 z + /3 k 4 


0, 


where z = x^o/lo — 1 > 0. Once again the solution is z — 1, or 


x = C 



Substituting into C(k), we have 


C(k) = 4 


7o/o 1 
Hq 1 + (3 k 4 ' 


(5.83) 


(5.84) 


5.3 Discussion 

As far as we are aware, the literature on statistical closures has neglected any kind 
of detailed examination of stability of the steady states. We believe the proof here of 
linear stability is the first such result obtained. In Appendix F, we provide a similar 
proof of stability for the Kraichnan-Spiegel closure (allowing for linear drive), which 
encompasses the Leith diffusion closure. A similar proof for EDQNM remains elusive, 
despite significant effort spent. (Orszag 197 ) clearly believes solutions to EDQNM to 
be stable. His arguments are compelling, although he was only considering turbulent 
drive due to external forcing, whereas we want to allow the case of arbitrary linear 
drive. It is not surprising that we have been unable to find a stability proof for 
EDQNM. EDQNM is far more complicated and may well allow for solutions which 
are linearly unstable in certain situations. A proof may not exist, or it may be beyond 
the limits of our imagination. 
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This closure for homogeneous turbulence could be extended with the appropriate 
terms for inhomogeneous interactions (which would mostly amount to pasting in 
terms from the CE2 equations). 
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Chapter 6 

Suggestions for Future Research 


Although this thesis has answered some questions, it has raised many new ones. Our 
theoretical analysis has taken place in the simplest possible setting, a 2D infinite 
(or periodic) system driven by white-noise external forcing. Naturally, one might 
wonder how to extend our analysis to more complicated, more realistic systems. For 
example, what happens in a realistic, physical geometry like a tokamak? If the 
system is driven by an intrinsic instability rather than external forcing, is there any 
qualitative difference? Do our previous results still hold in these instances? If not, 
why not, and how should the analysis be modified? 

We present a few of the issues in some detail, along with some ideas on how to 
proceed. Some of these issues could be studied within the CE2 formalism, while 
others would require more sophisticated statistical closures. This chapter necessarily 
includes some speculation in order to offer possible fruitful research directions. 

6.1 Using CE2 

There are several directions for future research even within the CE2 framework. First, 
one could perform some quantitative studies. For example, it has yet to be determined 
how the ZF length scale depends on other scales in the problem. For our numerical 
work we have done only two things. We have taken the deformation radius to be 
infinite, in which case the forcing length scale is the only external scale in the problem 
and sets the size of the ZFs. And we have taken the deformation radius to be of the 
same size as the forcing scale, in which case the ZF size must inevitably be similar to 
both. Parameter scans should be performed where the deformation radius and forcing 
scale are varied independently to determine their affect on the ZF size. Additionally, 
it would be interesting to study how large scale vs. small scale dissipation affects ZF 
saturation. For these studies, one might use DNS in addition to CE2. 

Second, one could attempt to understand in detail the problem with the Newton’s 
method used to solve for the ideal states numerically in Section 4.1. Near the in¬ 
stability threshold there was no issue, but far from threshold multiple solutions were 
appearing to the equations. The Newton’s method inevitably got stuck on nonphysi¬ 
cal solutions. Solving this problem would be worthwhile because our direct method of 
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solution of ideal states is otherwise limited to being near the threshold where the ZF 
is weak, and we cannot calculate the full stability diagram. There are a couple ways 
we envision proceeding. One might try using better numerical continuation methods 
( Vllgower and Georg 2003). These might do a better job of staying on the desired 
branch of physical, realizable solutions than the simple continuation method we have 
used. Additionally, others have used CE2 numerically with no problem (Farrell and 
Ioannou 2007, Tobias and Marston 2013); the difference between those methods and 
ours is that our method does not use a time evolution and excludes subharmonics. 
Therefore, one might try some kind of hybrid method involving both time evolution 
and Newton’s method to find the fixed point; the time-evolving method would help 
ensure a realizable solution. One could also include subharmonics in our calcula¬ 
tion. One other difference between our method and other numerical CE2 work is the 
use of an alternative coordinate system: we use the sum and difference coordinates 
y = y i — 2/2 and y = \{yi + 2 / 2 ) instead of y\ and y 2 - Using y 1 and y 2 has the advantage 
that the equations are in a form suitable for the Fast Fourier Transform, but the use 
of y and y allows us to change the ZF wavenumber q in tiny steps without changing 
the “box size” at the same time. 

The Hasegawa-Mima equation in periodic slab geometry omits a great deal of 
physics. We would like to understand zonal flows in toroidal devices such as tokamaks 
and stellarators. Some of the complications introduced are the magnetic geometry 
and linear instability. While linear instability is a topic we discuss in Section 6.2.2, 
the magnetic curvature leads to the existence of geodesic acoustic modes (GAMs). 
GAMs are modes with a zonally symmetric electric potential, like zonal flows, but are 
distinguished from zonal flows mainly in two ways: 1) GAMs oscillate at a frequency 
uj ~ c s /R, where c s is the acoustic speed and R is the major radius, and 2) GAMs 
are associated with a density perturbation that has sin 6 dependence, where 6 is 
the poloidal angle (Itoh et al. 2005, Winsor et al. 1968). Given how much we have 
learned about zonal flows using CE2, we have cautious optimism that something could 
be learned about GAMs as well. Besides for toroidal plasmas, linear plasma devices 
such as LAPD or CSDX may also provide a testbed and a window of understanding, 
especially for how shear flow interacts with turbulence. Linear devices are easier 
to analyze theoretically because of their simpler magnetic geometry and azimuthal 
symmetry. In linear devices, shear flow is often controlled through externally-applied 
potentials, although sometimes spontaneous shear flow emerges ( barter and Maggs 
2009, Holland et al. 2006, Tynan et al. 2006, Yan et al. 2010a,b, Zhou et al. 2012). 

I 11 the geophysical context, one obviously would want to know how these results 
extend to a rotating sphere. The /3 plane we have been using is merely an approxima¬ 
tion to the rotating sphere. We have been emphasizing the role of symmetry breaking, 
but moving to the surface of a rotating sphere destroys the north-south translational 
symmetries associated with a /3 plane. Do any of these results apply to zonal flows in 
spherical geometry? Although this question should be studied in detail, we offer one 
possibility. Due to the latitudinal variation of the Coriolis parameter, the turbulence 
is always inhomogeneous on the sphere. A transition from homogeneous to inhomo¬ 
geneous turbulence is not the right description, but perhaps some type of transition 
may still occur. Besides for the development of inhomogeneity, another aspect of the 
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bifurcation on a f3 plane is the spontaneous formation of a mean field, i.e., the zonal 
flow. We suggest that this mean-field generation may persist for flow on a rotating 
sphere, and would be observable as a control parameter is varied. The zonal flow 
still behaves as an order parameter in this more general type of scenario. This idea 
has some support, as numerical simulations appear to have observed this behavior 
as the rotation rate is increased from zero (Nozawa and Yoden 1997). Additionally, 
CE2 has been used to simulate turbulence on the rotating sphere, and ZFs have been 
observed within that framework ( vlarston et al. 2008, Tobias et al. 2011). Therefore, 
a future line of investigation could be to use CE2 to study zonostrophic instability on 
the sphere. This could be done numerically or possibly analytically by using equiv- 
ariant bifurcation theory (bifurcation theory for dynamical systems with symmetry) 
(Golubitsky et al. 1988). Qualitative insight could be gained into the structure of the 
unstable eigenfunction, including the direction of the equatorial jet. 

Finally, one could build upon the connection between zonostrophic instability and 
modulational instability described in Section 3.3 to improve our understanding of 
both. CE2 can be used to generalize modulational/secondary instability to more 
general background spectra. CE2 offers an alternative perspective into the physics 
of coherent-structure formation. It would be interesting to determine if CE2 can 
reproduce modulational/secondary instability when the eigenmodes are not Fourier 
modes (i.e., if nonperiodic boundary conditions are used). 

6.2 Other Statistical Formalisms and Closures 

In the theoretical study of turbulence, one line of approach is to examine statistically 
averaged quantities. That is the approach we have taken in this thesis, and it is dis¬ 
tinct from laboratory experiments or direct numerical simulation. In the statistical 
approach, one is interested often only in calculating second-order statistical quantities 
such as energy and transport, and so the closure problem arises for third-order terms. 
Many statistical closures of this type, which approximate the third-order terms in 
some way, have been studied in depth, including the Direct-Interaction Approxima¬ 
tion (DIA) and the Eddy-Damped Quasi-Normal Markovian closure ( 3owman and 
Krommes 1997, Bowman et al. 1993, Kraichnan 1959, Krommes 2002, Orszag 1977). 
These closures have several important properties. First, they conserve the same 
nonlinear invariants as the original dynamical equations through the same triadic 
mode-interaction structure. Second, they ensure statistical realizability. This means 
that statistical quantities are well-behaved under time evolution, so certain statistical 
constraints are guaranteed to be satisfied. For instance, realizability prevents energy 
from becoming negative. Some closures that do not respect realizability experience 
negative energies, an unacceptable flaw ( )gura 1962a, b). 

CE2, as previously discussed, can be categorized as a type of statistical closure. 
It is a closure for the one-time, two-point correlation function and allows for inhomo¬ 
geneous turbulence. CE2 is particularly simple, since the closure technique involves 
nothing more than neglecting the unknown terms. This means that CE2 totally 
ignores eddy self-nonlinearities, which are responsible for the traditional cascades. 
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Using one-time correlation functions rather than the more general two-time functions 
means that CE2 also lacks certain time-history information and loses some of the ef¬ 
fects of wave propagation (Krommes and Smith 1987). To incorporate these physical 
effects, as well as to achieve greater quantitative accuracy, the effect of eddy self- 
nonlincarities and time-history information must be retained in some way through a 
more sophisticated closure like those described above. 

Another type of approach does not focus solely on second-order or nth-order sta¬ 
tistical quantities, but uses the full probability density function (or functional). This 
is the approach taken by Bouchct et al. (201 ), who used it to rigorously justify the 
quasilinear approximation in the barotropic vorticity equation in the large 7 limit. 

The averaging procedure to obtain the CE2 equations from the QL equations mer¬ 
its further discussion ( arker and Krommes 201 ). We used a zonal average, but for 
CE2 or other formalisms, other types of averages may be used. Under appropriate as¬ 
sumptions, which always include some kind of ergodicity assumption, multiple choices 
of average will lead to the same final equations. For instance, zonal (Srinivasan and 
Young 2012), short-time (Bakas and Ioannou 2011), and coarse-graining (Bakas and 
Ioannou 2013a) averages have been discussed. The ergodicity assumption allows one 
to transform the average over the random forcing into a deterministic quantity. One 
can also discuss things in terms of an ensemble average, in which case an assumption 
of statistical homogeneity in the zonal (a:) direction is made, but inhomogeneity is 
allowed in the nonzonal ( y ) direction. In this case, ergodicity is not required in order 
to derive the CE2 equations, but it becomes necessary if one wants to interpret the 
solutions of the equations as having anything to do with the behavior of an individual 
realization. 

When using the ensemble average, Kraichnan pointed out in the context of ther¬ 
mal convection that the definition of the statistical ensemble is somewhat subtle for 
the situation of spontaneous symmetry breaking ( kraichnan 1964 ). Because of the 
translational symmetry, the zonal jets have no preferred location and are presumably 
equally likely to form with any particular phase. One choice of the statistical ensem¬ 
ble encompasses all possible realizations consistent with the prescribed parameters, in 
which case the ensemble itself is statistically homogeneous and any ensemble-averaged 
quantity must be homogeneous also. Therefore the average yields zero mean ZF (and 
then the ZF must be described as a fluctuation), despite the fact that each individual 
realization has a nonzero ZF. This was the procedure followed in Krommes and Kim 
(2000). Another possibility is that the ensemble might consist only of the realizations 
for which the zonal jets have a particular phase. The latter interpretation is the one 
that yields the CE2 equations identical to those obtained by zonal averaging. With 
the former ensemble, the ergodic assumption is invalid, since an individual realization 
is no longer mixing throughout the full set of realizations of this ensemble. This is 
consistent with the fact that the ensemble-averaged behavior is not equivalent to the 
behavior of an individual realization. 
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6.2.1 Development of Systematic Closures for Inhomoge¬ 
neous Turbulence 

Historically, the majority of analytical theories of statistical turbulence assume homo¬ 
geneous statistics, where the statistics of turbulent quantities do not depend on po¬ 
sition. Relatively little effort has been devoted to inhomogeneous statistics. Progress 
developing inhomogeneous closures has been limited and is one area for future re¬ 
search. 

One proposed way to go beyond CE2 is to use third-order cumulants in a CE3 
framework, where fourth-order cumulants are neglected ( fbbias and Marston 2013). 
That could be useful when eddy-eddy nonlinearities are a small perturbation. But 
this approach has problems because CE3, unlike CE2, is not realizable; it must be 
patched up in an ad-hoc manner. 

A few systematic inhomogeneous closures exist, mostly stemming from Kraichnan. 
One is the full, inhomogeneous DIA ( kraichnan 1964 ). Kraichnan also proposed a 
simpler DIA variant called the diagonalizing DIA (Kraichnan 1964a). More recently, 
the diagonalizing DIA has been generalized into the quasi-diagonal DIA (Frederiksen 
1999, O’Kane and Frederiksen 2004), but these “diagonal” DIA closures approximate 
the interaction between the mean held and the fluctuation. That approximation would 
affect the stability properties of the ZF in ways currently unknown. Additionally, 
an inhomogeneous Markovianized closure exists in the test-held model ( kraichnan 
1972), but it is not statistically realizable in the presence of waves ( Bowman and 
Ivronnnes 1997, Bowman et al. 1993). A homogeneous realizable test-held model 
exists (Bowman and Kronnnes 199" ), but as of yet there is no version that is both 
realizable and inhomogeneous. More work along these lines needs to be done. 

6.2.2 Systems with Intrinsic Instability 

The statistical closures described above are intended to more faithfully represent the 
eddy-eddy nonlinearities than CE2 does. This can be important for more than mere 
quantitative accuracy. We can imagine at least one situation for which it is crucial 
to retain the eddy self-nonlinearities: a system with linear instability. Linear in¬ 
stabilities in plasmas are common, such as the ion-temperature-gradient instability. 
And the oceans are baroclinically unstable. Numerical simulations of the Modihed 
Hasegawa-Wakatani system have clearly demonstrated the symmetry-breaking bifur¬ 
cation of ZF generation from homogeneous to inhomogeneous turbulence (Numata 
et al. 2007). In order to describe this transition, a model must allow for an equi¬ 
librium of homogeneous turbulence. But in a quasilinear (QL) or CE2 description, 
if no ZFs are present then there are no nonlinear interactions, and it is impossible 
for a linear instability to saturate. With linear instability present, a QL description 
has no homogeneous equilibrium. Retaining the eddy self-nonlinearities is required 
to allow a statistical equilibrium of homogeneous turbulence, which can then undergo 
zonostrophic instability to generate zonal flows. 

A schematic of the different possible regimes as a function of parameter space is 
sketched in Figure 6.1. The transition from region 1 to region 2 gives the transition 
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Figure 6.1: Hypothetical schematic of three possible parameter regimes in an ion 
temperature gradient (ITG) system, as a function of two parameters, the zonal flow 
damping rate fj, z and the temperature gradient R/Lt■ The point A indicates the 
codimension-2 bifurcation point where regions 1, 2, and 3 intersect. 


to homogeneous turbulence as linear instability becomes active. The transition from 
region 2 to region 3 is the zonostrophic bifurcation studied in great detail in Chapter 3, 
where homogeneous turbulence becomes inhomogeneous as zonal flows are born. The 
transition from region 1 to region 3 is not understood at this point. The point A 
indicates the coclimension-2 bifurcation point where regions 1, 2, and 3 intersect. A 
bifurcation analysis about the point A might be interesting. 

However, this sketch may be too simplistic for even the least complicated plasma 
turbulence systems studied. For one, we have assumed there is no subcritical turbu¬ 
lence 1 . We have also let the dissipation parameter of the zonal flows, p 2 , be controlled 
independently from other parameters. But in the Modified Hasegawa-Wakatani sys¬ 
tem, the dissipation is not so simple, and it may not be controlled directly. Instead, 
much of the dissipation arises from the coupling of the zonal flows to nonzonal modes, 
which then suffer from resistive damping (Hatch et al. 2011a,b, Makwana et al. 2012, 
2011, Terry et al. 2006). The dissipation is determined nonlinearly after saturation 
by all the mode couplings. That kind of scenario will have to be studied in detail. 
The sketch we offered in Figure 6.1 is just a beginning. And there may be other types 
of regimes and transitions that we are yet unaware of. 

1 Subcritical turbulence refers to turbulence that is sustained even when the base state is linearly 
stable. 
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The Dimits shift is another aspect of certain linearly unstable magnetically con¬ 
fined plasmas ( limits et al. 2000). The Dimits shift has been a phenomenon of high 
interest ever since it was discovered numerically, and there is as yet no experimental 
evidence for it. In the Dimits shift, turbulence and turbulent heat transport are sup¬ 
pressed even beyond the linear marginal stability boundary. In other words, when 
the ion-temperature gradient was increased to just beyond the critical value for linear 
instability, no turbulence was observed. This unexpected behavior was attributed 
to the suppression of turbulence by nonlinearly-generated zonal flows. As the ion- 
temperature gradient was increased even further, eventually turbulence and turbulent 
transport would develop (possibly because the zonal flows suffer their own instability 
and can no longer effectively suppress the turbulence). This upshift from the linear 
stability boundary to some other boundary is termed the Dimits shift. 

With collisionless ZFs, a Galerkin-truncated ITG system of just 10 modes was 
found to exhibit a Dimits shift ( folesnikov and Krommes 2005a, b). However, it is 
unclear what exactly can be learned from that calculation, because the behavior of 
the system was sensitive to the number of modes retained in the truncation. Perhaps 
an analysis that retains the full spatial dependence, through the inhomogeneous sta¬ 
tistical closures we have been describing, would lead to more regular and well-behaved 
behavior and improved understanding of the Dimits shift. One possibility, suggested 
by the pattern formation framework, is that in the Dimits shift regime, steady zonal 
flows exist within some stability balloon. But at large enough profile gradients, any 
steady zonal flow becomes unstable, leading to rapidly fluctuating zonal flows and 
reduced suppression of turbulence. This scenario would be consistent with the ideas 
of fogers et al ( !000). 

The Dimits shift is not understood theoretically. Many studies of it use collision¬ 
less ZFs, but not all ( likkelsen and Dorland 2008). What can be said is that the 
Dimits shift involves a transition that includes the generation of zonal flows. This is 
a type of behavior similar to the zonostrophic bifurcation that has been successfully 
described by CE2. It is possible that the statistical framework with inhomogeneous 
turbulence may be similarly successful in describing the Dimits shift. To follow this 
route, one would want to find the simplest system that exhibited a Dimits-shift-like 
behavior. For example, does the Dimits shift require the effects present in gyroki- 
netics, or can it be adequately captured in a fluid description? Is toroidal magnetic 
geometry essential, or is there a simpler geometry that possesses sufficiently similar 
behavior? A minimal model would make the analysis and physics as transparent as 
possible. 

The zero-dimensional phenomenological bifurcation model of zonostrophic insta¬ 
bility, presented in Section 3.1, can be modified for the case of a linear instability by 
the addition of terms representing the eddy self-nonlinearity. For instance, following 
the example of a nonlinear closure with quadratic terms, one might have 

W h = 7 W h - ixWl + FWl - aWiZ, 
m = 7 Wi - fiWhW, + v W h z , 
z = —vz + aWi, 
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( 6 . 1 a) 

(6.1b) 

( 6 . 1 c) 


where we have assumed that the incoherent forcing term F does not appear in the 
inhomogeneous equation. (This model has the problem that in some circumstances 
W h can become negative and blow up.) The linear instability term 7 might be related 
to the temperature gradient R/L T . We can absorb F into /i and thus write 

W h = 7 W h - HhWl - cvWiZ, (6.2a) 

Wi = 7 Wi - HiWhWi + 7}W h z, ( 6 . 2 b) 

z = — p, z z + aWi, ( 6 . 2 c) 


where /i, > /x / 1 . The zero state is linearly unstable if 7 > 0. This model has a 
homogeneous equilibrium at Wh = 7 / Hh- Its linear stability can be checked easily. 
The condition for zonostrophic instability is 


arj fa 
HhHz h/i 


(6.3) 


This simple model has a structure similar to that in Figure 6.1. When 7 > 0, the zero 
state is unstable, with a zonostrophic boundary depending on the value of Hz (with 
all other parameters fixed). In this simple model, the zonostrophic boundary has no 
dependence on 7 . The model is also not complicated enough to have a Dirnits shift. 

The closure for homogeneous statistics presented in Chapter 5 could be extended 
to inhomogeneous statistics as well. That closure has a nonlinear damping term and 
so can handle intrinsic linear instabilities. Since it is rather simple, analytic progress 
might even be possible, e.g., in a bifurcation analysis. 


6.3 Other Gaps in Knowledge 

In pattern-forming systems, the simplest theoretical approach is to eliminate bound¬ 
aries and use an infinite or periodic system. That was the initial approach taken in 
Rayleigh-Benard convection, and that is the approach taken here. However, bound¬ 
aries are actually quite important. For instance, one might expect that if a system 
gets very large, then far from the boundaries, the boundaries have little effect. But in 
the amplitude equation (3.60), a prototypical pattern-inhibiting boundary condition 
A — 0 has a profound effect on the possible wavenumbers of the pattern even far 
from the boundary. Instead of a band of stable, stationary solutions as in the case 
of infinite boundaries, an A = 0 boundary condition in a semi-infinite system forces 
the pattern wavenumber to be unique and equal to the critical wavenumber (Cross 
and Greenside 2009). Some boundaries can suppress the amplitude of patterns, and 
others can enhance pattern formation. 

In toroidal and cylindrical plasma devices, boundaries exist and have a major in¬ 
fluence on the plasma’s behavior. In toroidal plasmas, the magnetic geometry plays a 
dominant role in determining the character of the turbulence, and we should expect 
that the magnetic geometry and especially the separatrix influence the generation 
and characteristic of zonal flows. Systematically understanding the geometry and 
boundary effects is a major open area for study and will rely heavily on simulations. 
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Some initial work has been done in terms of examining how various stellarator con¬ 
figurations affect ZFs (Xanthopoulos et al. 2011), but far more work needs to be 
done. 

In many simulations in toroidal geometry, ZFs are observed to be non-steady. This 
fluctuating behavior is distinct from the steady ZFs we have been assuming in the 
theoretical analysis in this thesis. If the time scale of the ZFs 7 fluctuations are long 
compared to that of the turbulence, then perhaps an assumption of steady ZFs is an 
acceptable lowest-order approach. But ZFs have sometimes been seen to fluctuate on 
the same time scale as the turbulence, in which case the theory developed here is not 
directly applicable. 

How can we use any of this knowledge to benefit experiments, or even to talk in 
a language that experimentalists understand? In the geophysical context, possibly. 
Given the numerous discoveries of exoplanets and the ever-more sophisticated obser¬ 
vational methods, we someday might encounter an exoplanet gas giant that has no 
zonal jets. This would contrast with the gas giants within our solar system, which 
all have zonal jets. A fundamental theoretical understanding of the zonostrophic 
bifurcation is key to puzzling out how various factors impact zonation. 

Plasmas, on the other hand, are so messy and complex that we currently see no 
direct way for the theory to be directly compared with experiment. The Hasegawa- 
Mima equation neglects many, many physical effects. We discovered some general 
principles in the 2D slab geometry, but it is unclear if those survive in toroidal geom¬ 
etry. Even the cylindrical plasma devices, with their simpler magnetic geometry, are 
so small that radial boundary conditions inevitably have a strong influence. 

To us, the way to proceed to develop this theory for usefulness to plasma physicists 
is twofold. One direction is to increment in complexity, step by step. For instance, 
eddy-eddy nonlinearities can be added to handle linear instabilities. The theory 
should be constructed in cylindrical geometry, then in toroidal geometry. GAMs 
should be investigated. Kinetic effects might be added. A worthy goal would be 
to try to identify and understand the Dimits shift in a simple model, as explained 
above. The second direction goes hand-in-hand with the first, and that is to firm 
up the theory with numerical simulation. We believe that many of the principles 
that we have found from the QL approximation to the Hasegawa-Mima equation will 
apply in many other cases. It appears generic that steady zonal flows are generated 
in slab geometry in plasmas. Some gyrokinetic ITG simulations in slab geometry 
have seen steady zonal flows ( tatch), in which case the pattern formation principles 
ought to apply. Detailed comparisons of such simulations with theoretical predictions 
will undoubtedly lead to progress. It is only by laying the groundwork that we as 
a community will be able to construct the elaborate theoretical towers required to 
understand plasma turbulence. 

Finally, one area of high interest, which was originally to be one of the questions 
considered in this thesis but was barely touched on, is how zonal flow suppresses tur¬ 
bulence. Multiple explanations have been given, but there is no firm theoretical basis 
for which to understand and compare them. Since the pattern formation approach is 
new in the field of zonal flows, it provides a novel way to attack this problem. 
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Appendix A 

Derivation of CE2 in Real Space 


Here we provide the details of the derivation of the CE2 equations (2.21). This 
procedure follows that by Srinivasan and Young ( 012). 

We begin from the QL system (2.17), which we restate here: 

d t w' + {[/V 2 + p- [(<% - L- 2 )U]}d x iP' = £ - /W - v{-l) h V 2h w’, (A. la) 

[dt + y + u(-l) h d 2h ] (l - a ZF L d 2 d- 2 )U{y) + dyV^J y = 0, (A.lb) 

The covariance of the white-noise forcing £ is defined to be 

(£(ar, 2 / 1 , ti)f (a; 2,1/2, t 2 )) = F(x i - x 2 , yi - y2)S(h - 1 2 ) (A.2) 

The forcing is taken to be homogeneous in space such that its statistics only depend 
on the spatial difference x x — x 2 = (x\ — x 2 ,yi — y 2 ). 

Define 

W(x 1 ,y 1 ,x 2 ,y 2 ) =w\x 1 ,y 1 )w , (x 2 ,y 2 ) (A.3) 

(taken at the same time t). Averaging W over x\ holding x 2 fixed (or vice versa) 
gives zero, by definition. Instead we define 

W(x 1 ,y 1 ,x 2 ,y 2 ) = W (x x - x 2 l y x - y 2 \ \{x x + x 2 )) = W{x, y { , y 2 \ x), (A.4) 

with the sum coordinate x — \{x\ + x 2 ) and the difference coordinate x = X\ — x 2 . 
At a later point, we will also switch to sum and difference coordinates for y. Now, 
define 

1 f Lx _ — 

W(x,y 1 : y 2 ) = — dx\ x W(x,y 1 ,y2\x) 

Jo 

1 f Lx _ 

= — / dx\ x w'(x 1 ,y 1 )w'(x 2 ,y 2 ), (A.5) 

X'x Jo 

where L x is some averaging length. This averages the product w'(xi,yi)w'(x 2 ,y 2 ) 
holding the separation X\ — x 2 fixed. This zonal average presumably smooths rapidly 
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fluctuating quantities (in space and time). Similarly, we can define 


^(x,yi,y 2 | x) = i//(xi,yx)‘il/(x 2 ,y 2 ), (A.6) 

'&(x,y 1 ,y 2 ) = — / dx\ x 'S>(x,y 1 ,y 2 \ x) = — / dx\ x ^'(xi 1 y 1 )'i/j'(x 2 ,y 2 ). (A.7) 

X>x Jo X>x Jo 

We can relate W and T. Recall that w'(x,y ) = V 2, ip'(x,y) = (V 2 — L^ 2 )jp'(x,y). 

Then 

w'ix^yx) = [d 2 xl \ X2 + d 2 yi — Lj 2 ]i/j'( xi, yi), (A.8) 

w'(x 2 ,y 2 ) = [d 2 X2 \ Xl + d 2 m - L~ 2 ]'i/j'(x 2 ,y 2 ). (A.9) 

Now, use d Xx = d x + \d X) and d X2 = —d x + \dx. Substituting these relations into 
(A. 3), we have 

W(x, y lx y 2 ) = j- [ dx\ x [d 2 x + d xx + \&l + d 2 yi - L~ 2 ] 

J 0 

o [d 2 x - d xx + \dl + d 2 yi - L~ 2 ] ^\x u yi)^\x 2 , y 2 ). (A.10) 

By the assumed periodicity in x (or other assumption), the dx terms vanish. Dehne 



(A.ll) 


(A.12) 

for j = 1,2. Then, we see that 


W(x,y u y 2 ) = V?V 2 ^(x,J/i,S/ 2 )- 

(A.13) 

In shorthand notation, we also write 


W(x,yi,y 2 ) = w[w' 2 , 

(A.14) 


where w' = w'(xj , yj ) and the overbar means spatial average holding x = x\— x 2 fixed. 
Similarly, for the velocity correlation tensor, one finds (with u = v x and v = v y ) 


V ij (x 1 ,x 2 ,y) 


u'iU 2 u\v£ \ 

U 2 v[ v[v' 2 J 


3 3 

U yi U y2 

—3 3 

w l)2 


d x d yi \ 

-dl) 


^{x,yi,y 2 ). 


(A.15) 


Because the choice of denoting one point as xi and the other as X 2 is arbitrary, all 
correlation functions have the exchange symmetry (Srinivasan and Young 2012) 


W(x,y 1 ,y 2 ) = W(-x,y 2 ,yi). 


(A.16) 


Now, we derive an evolution equation for W. From (A. 14) we have 


d t W = (d t w[)w2 + w[(d t W2). 
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(A.17) 









Substituting in from (A. la) and applying the averaging, one eventually finds 


d f W + (v 2 2 L ± - V?L 2 ) d x V = -2 /iW - v{-l) h (Vf + Vf) W 


^W' 2 + W'^2, 


(A.18) 


where 


Lj = UjV* + ((3 — u'j), 
Uj = Ufa), 
u" = d 2 yj U( yj ), 

for j = 1, 2. For later use, notice that 


(A.19) 
(A.20) 
(A.21) 


=2, —2 . 


^rlU=2 


. ^2 


V 2 L 1 ~V 1 L 2 )d x * = (U 1 ~U 2 )d x W + (/3^C7 1 )V 2 -(/3-C/ 2 )V 1 d x *. (A.22) 


Now, we switch to using sum and difference coordinates in y, with y = r/i — y 2 
and y = (y 1 + y 2 )/ 2, with d y] = d y + ^<9^ and d y2 = —d y + We write 

W 7 0, 2 / 1 , 1 / 2 ) = VF(x, y | j/). (A.23) 

In terms of y and y, the Laplacians are 

v? = V 2 + d,d s + 1S|, (A.24) 

V 2 = V 2 - dydy + (A.25) 

where now V 2 = <9 2 + d y is the “separation” Laplacian. We also define 

V 2 = V 2 - L~ 2 . (A.26) 

The symbols V 2 and V" were used in slightly different context in the fluctuating 
amplitude equations, but now we reuse them purely in the averaged equations and 
the meaning should be clear. From (A. 13), we can relate W and T in the new 
coordinates, 

W(X, v\y)= (V + dA+ \ai\{v 2 - d y d,+ l -^(x,y\y). (A.27) 

Using the sum and difference coordinates, the evolution equation for W, (A.18), 
becomes after some algebra 

d t W(x, y\y) + {U + - U.)d x W - (i U" + - ff)(V 2 + \d 2 )d x ^> 

- [2/5 - (u" + + u"_ )]dydyd x ^ = (6w 2 + w[&) x - 2 fiW - 2 uD h W, (A. 28 ) 
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where now U± — U(y ± \y), U± = U"(y ± \y) — &zFL d 2 U(y ± \y), and Dh is the 
hyperdiffusion operator, given by 


Dh = t- 1 )^ 


1 h 


+ (dy + \®y) + ^ + (< 9 ; 



(A.29) 


We must now compute the term resulting from the external stochastic forcing, 
i\w' 2 + w[£ 2 - We make an ergodic assumption such that a zonal average is equivalent 
to an ensemble average over the realizations of the forcing, 


6^2 + w i& = (6^2 + w i&)- (A.30) 

With this assumption, the desired term can be calculated exactly. The assumption 
that the forcing is white noise (delta-correlated in time) is also crucial. From (A. la), 
we can write 

w\x 2 ,y 2 ,t) = w 0 (x 2 ,y 2 ,t 0 ) + f dt'N(t')+ f dt'£(x 2 ,y 2 ,t'), (A.31) 

Jt 0 Jt 0 

where to < t and N contains all the appropriate terms. The ensemble average 
(£(xi, 2 / 1 , t)w'(x 2 , y 2 , t)) becomes 

wo{x 2 ,y 2 ,t 0 ) + [ dt'N(lf)\£(x 1 ,y 1 ,t)\+ f dt'(^(x 1 ,y 1 ,t)^(x 2 ,y 2 ,t 1 )). (A.32) 
Jt 0 J / Jt 0 

The hrst average vanishes because the fields wq and N at times prior to t are un¬ 
correlated with the random forcing at time t. The second average is given by the 
definition of the forcing (A.2). One is left with the integral 

(£ 1 ^ 2 ) = F(xi ~ x 2 ,yi - y 2 ) f dt'5(t-t'). (A.33) 

Jto 

The integral over the delta function is somewhat subtle because t' — t occurs exactly 
at the endpoint, but it gives exactly This can be seen intuitively because any phys¬ 
ical correlation function must be nonsingular and symmetric about its time argument. 
Thus half of the ‘weight’ of the correlation function sits at t < t! and the other half at 
t > t 1 . If one considers white noise as the limit of some process with finite correlation 
time, then one must conclude that only half of the correlation function is integrated 
over, leading to the value of the integral as Similarly, it is not hard to check that 
(^£ 2 ) evaluates to the same result of ^F. 

Thus the evolution equation for W becomes, finally, 

dtW{x, y\y) + (U+- U.)d x W - {U" + - u"_) (V + <9 X T 

- [2P-(UI + U^)\dyd y d x <!! = F(x,y)-2fj.W -2uD h W. (A.34) 
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The Reynolds stress term in the equation for the zonal flow can be written in 
terms of the eddy correlation function. In sum and difference coordinates, mean- 
square quantities are obtained by evaluating correlation functions at zero separation, 
i.e., by setting (x,y) = 0. For example 

w{x 1 ,7/iV(zi, 2/1 ) = W (x, 2 / 1 , 2 / 1 ) = w(0, 0 I 2 / 1 ) (A.35) 

(with y = yi). From (A. 15), we have 

u[v 2 + = 2<9 x <9 y \F (A.36) 

Evaluating at x 2 = aq and 2/2 = 2/1 = 17, so that £ = 0, y = 0, we have 

Wv'(y) = d x d y ^(0 , 0 | y). (A.37) 

Thus, as a function of y, the mean flow equation (A. lb) can be written as 

[3* + y + v{-l) h df ]iu(y) + dyd x d y ^(0,0 I y) = 0, (A.38) 

where 1 = 1 — azpLy 2 d^ 2 . 

Equations (A.34) and (A.38) form a closed system called CE2 (second-order cu- 
mulant), along with (A.27) relating W(x,y \ y) and \I7( x,y \ y). 
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Appendix B 

Correlation Function 
Corresponding to a Wave 


We consider in this section the one-time, two-point correlation function corresponding 
to a wave. First we consider the general case of a superposition of waves. Let 

ip'(x,y,t) = 2^2c k cos(k x x + k y y - uj k t + 0 k ). (B.l) 

k 

Then, letting xp( = ip'(xi,y\,t) and 'iJ/ 2 — ' l P'(x 2 ,y 2 ,t), we have 
0 i 02 = 2c k c k'{cos [\{k x + k' x )x + (/c x - k' x )x + \(k y + k' y )y + (, k y - - % k ,] 

k, k' 

+ cosfK^ - k' x )x + (A:* + A£)x + \{k y - k' y )y + + k' y )y - z^,\ }, (B.2) 

where x = x\ — X 2 , x = |(xi + X 2 ), and z kk , = (c^k ± Wk >)t — (0k ± 0k')- Using a zonal 
average, the correlation function is obtained by integrating over x with x held fixed: 

^{x,y\y) = j-( dx\ x xp[xp 2- (B.3) 

X'x Jo 

The first cosine vanishes unless k' x = k x , while the second cosine vanishes unless 
k' x = —k x . For simplicity assume all the kxik x > 0. Then we are left with 

^(x,y | y) = y t 2c k c k ' cos[k x x+l(k y +k' v )y+(k y -k y )y-(uj k -u k >)t+4> k -(l) k '\. (B.4) 

k.fc' 

If we separate out in the sum the terms for which k' y = k y , then we have 

V(x,y | y) = E 2 c k cos (k x x + k y y) + EE 2 ckCk' cos [k x x 

k k k' y jXiky 

+ \{k y + k' y )y + (k y - ky)y - (ca k - uj k >)t + 0 k - 0k'] • (B.5) 
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It can be verified by substitution that this is a solution to the unforced, undamped 

—2 

CE2 equations without zonal flow, d t W = 2 /3dyd y d x ^ (and using Wk = —k x /3/k ). 
We see that the first term of (B.5), which corresponds to the covariance of individual 
waves, is unchanging in time and homogeneous in space. But in the second term, 
waves with different k y give rise to a correlation function that oscillates in time and 
has y dependence. This is a manifestation of the coherent beating between waves. 
There is no decorrelation mechanism present; that requires nonlinear physics. 

One can imagine using another averaging procedure instead of the zonal average. 
With the zonal average, the only coherent structures allowed are zonally symmetric. 
One might also want to investigate zonally asymmetric structures, which precludes 
the use of a zonal average (Bakas and Ioannou 2013a). To study these more general 
coherent structures, the correlation function can be defined using a coarse graining 
in time or space (this approach typically requires the mean field and fluctuations to 
obey a scale-separation assumption) or an ensemble average. 

To illustrate an alternate derivation for a single wave, let 

V>'(x) = i/j 0 (e ip ' x-iwt + e~ ipx+iut ) . (B. 6 ) 


Then 


Vh"02 = i’l {e 2 ip ' x e~ 2iujt + e ipx + e~ ipx + e ~ 2 ip ' x e 2iut ) . (B.7) 

At this point, a coarse graining in time over an intermediate time between cu _1 and 
the timescale of the coherent structure eliminates the oscillating terms. Equivalently, 
one could perform a coarse graining in space over an intermediate scale between p ~ 1 
and the size of the coherent structure. Then, one obtains 

T = (e ip x + e~ ip x ) . (B.8) 

This T is homogeneous (independent of x). Its Fourier transform is 

'M k) = (27t) 2 V>o [<S(k - p) + 5(k + p)] . (B.9) 

The inclusion of the mode at — p as well as the mode at p is essential and arises from 
the reality condition. 
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Appendix C 

Derivation of the Amplitude 
Equation 


Here we derive the amplitude equation (3.60) directly from the CE2 equations (2.21) 
and verify the results numerically. First, we review the procedure for the perturbation 
expansion (Cross and Greenside 2009). Then we fill in the algebraic details. 

C.l Review of the Perturbation Expansion 

We limit ourselves to quadratic nonlinearity. Let (j) be an abstract vector, A be a 
linear operator, A be a bilinear operator, and F be external forcing. Any of A, N, 
and F may depend explicitly on the small parameter e. The basic equation is taken 
to be 

0 = A0 + N{(j>, <j>) + F. (C.l) 

Without loss of generality, N can be assumed to be symmetric in its arguments (if 
it is not, a new symmetrized operator can be defined and used instead). Given a 
nonzero equilibrium </> e , we change variables by letting (f> = (f> e + u to give 

0 = Lu + N(u,u), (C.2) 


where Lu = A u + 2 N(4> e , u). 

We take as given that at e = 0, the equilibrium </> e transitions from stable to 
unstable due to a perturbation with wavenumber q c . This calculation is motivated 
by the discovery of the zonostrophic instability, described in Section 3.2. Figure C.l 
depicts the schematic of the bifurcation. 

In performing the perturbation procedure, we use a multiple-scale expansion with 
slowly varying space and time scales. This is accomplished by introducing the slow 
scales Y = e l ^ 2 y and T = et, then letting dy —* dy + F^dy and d t —> d t + edy- Using 
these, we expand L = L 0 + F^Lx + eL 2 + e 3 / 2 L 3 + • • •, N = N 0 + e 1//2 Ah + ■ • •, 
and u = F^ 2 U\ + eu 2 + ■ ■ ■ ■ Expansion in F^ 2 rather than in e arises due to generic 
behavior of supercritical bifurcations. Collecting terms of the same order, we obtain 
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Figure C.l: Schematic of bifurcation. Top: The homogeneous equilibrium xh is stable 
(solid) for e < 0 and zonostrophically unstable (dashed) for e > 0. At e = 0, a new 
set of inhomogeneous equilibria xj appears; some of these equilibria are stable and 
some are unstable. Bottom: Growth rate A of perturbations about the homogeneous 
equilibrium xh as a function of the ZF wavenumber q. 




(C.3) 

(C.4) 

(C.5) 


the equations at 0 (e 1//2 ), 0 (e), and 0 (e 3 / 2 ): 

0(e 1/2 ) : 0 = L 0 ui, 

0(e): 0 = L 0 u 2 + Liu-! + N 0 (ui,ui), 

0(e 3/2 ) : 0 = L 0 u 3 + Liu 2 + L 2 ui + 2N 0 (ui,u 2 ) + Ni(ui,ui). 

At 0(e 1//2 ), the condition L 0 ui = 0 states that u\ is an eigenvector with a zero 
eigenvalue. Then u\ can be a linear combination of null eigenvectors with a to-be- 
determined amplitude. The reality condition on u restricts the form to be 

u x = A(Y, T)r + A(Y , T)*r *, (C. 6 ) 

where r ~ e iqcq (and its complex conjugate) are the right null eigenvectors. These 
eigenvectors are periodic in y with critical wavenumber q c , which is the hrst wavenum¬ 
ber to go unstable as e crosses zero. Given an inner product (•,•), then associated with 
the right null eigenvector is a left null eigenvector l of L 0 , such that (l, L 0 u) = 0 for 
any u. The y dependence of l will also be e* 9c U The amplitude A will be determined 
by nonlinearities occurring at higher order. 

At 0(e), we hrst note that LiU\ = 0 automatically. This is because q c is marginally 
stable at the instability threshold: given a dispersion relation A (q,e) as a function of 
wavenumber q and control parameter e, then both A(g c ,0) = 0 and d\/dq(q c ,0) = 0 
(see Figure C.l). The former equality yields L 0 Ui = 0 and the latter equality yields 
the condition L\U\ = 0. In order to ensure that a solution for 112 exists, a solvability 
condition obtained by taking the inner product with the left null eigenvector must 
be satisfied. This solvability condition is (/,L 0 m 2 + At 0 (u 1 ,u 1 )) = (l, N 0 (ui,ui)) = 0. 
Because l ~ e iqcV and N 0 (ui,ui) ~ 1 or e ±2tqcq due to the quadratic nonlinearity, this 
solvability condition is always satisfied. Thus, given that a solution exists, one may 
write U 2 as a linear combination of homogeneous and particular solutions: 

U2 = U2h + U2pi (C.7) 

where 

u 2h = A 2 (Y, T)r + A 2 (Y, T)*r*, (C.8) 

L 0 u 2p = -JVo(iti,ui). (C.9) 

Since we have not yet determined A, we must proceed to higher order. Another 
unknown parameter A 2 has been introduced, but we will not need it in order to solve 
for A. 

At 0(e 3 / 2 ), note that L 1 u 2 h = 0 for the same reason that LiU\ = 0 . Upon writing 
the solvability condition from (C.5), one finds that several terms vanish, leaving 

0 = (li, L 2 Ui) + (li, 2N 0 (ui, U 2 P )) • (C.10) 
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This is the desired partial differential equation which determines the amplitude A. 
Note that one never explicitly needs L\ or N 1 . 

C.2 Details 

We now apply this procedure to (2.21). For simplicity, we set the viscosity to zero, take 
infinite deformation radius, and cross the instability threshold by varying the strength 
of the forcing (rather than by varying the friction as in the main text); modification 
for other scenarios is obvious. Let the forcing be given by F(x,y ) = (1 + e)F 0 (x,y), 
where instability threshold is at e = 0. We shift variables relative to the equilibrium 
at Weqb = (1 + e)F 0 /2y } U eq b = 0. Explicitly, the abstract vector u consists of two 
components, u = {W, U}. Then we have the basic structure of (C.2) with 

[L{u)] w = -d t W(x,y,y) + 2/3d y d x dyD~ 1 W(x,y,y) - 2yW(x,y,y) 

- + \y) ~ u (y~ ly)) d xF 0 (x,y) 

-ai[U(y+ \y)-U(y-\y)}F Fa ( x ,y)Y (C.lla) 

[£(«)]u = -a,U(y) - d„a x d,D- 1 W(0 ,0,5) - yU(y), (C.llb) 

where W in (C.llb) should be evaluated at x = 0 and y = 0 after performing the 
derivatives, [-]w an d [•][/ refer to the W and U components of the abstract vector, 
and 

o = V 4 + + Laj. (C.12) 

Note that D commutes with d x , d y , and d y . 

For the nonlinear operator N(v, z ), with v = {vw, vu} and z = {zw, zu}, we have 
the unsymmetrized version lV uri , 

[N un (v,z)] w = —[vjj(y + \y) - v v (y - \y)]d x z w 

+ [d^v v (y + | y) - d*vu(y - ||/)](V 2 + \c^)d x D~ 1 z w 
- [d^vu(y + \y) + d^v v (y - \y)]d y d x d y D^ l z w , (C.13a) 

[A^ un (u, z)\ v = 0. (C.13b) 

The symmetrized operator is then given by 

N(v,z) = 1[JV“(b, 2 ) + JV'“( 2 ,b)]. (C.14) 

We now introduce the slow space and time scales. One subtlety that was not 
mentioned in the general procedure described above is that the Ui may need to be 
expanded in e. This occurs for two reasons. First, because the U(y ± \y) terms lead 
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to 

A{e l l\y±\y))=A{Y±\t 1 / 2 y)=A{Y)+ l -ye l / 2 d Y A+ l -y\d 2 Y A-V--- . (C.15) 

Second, the right null eigenvector r itself contains the differential operator dy (i.e., 
it depends on q ), which must be expanded in the multiple-scale procedure. It is 
extremely convenient to introduce these expansions at the outset so as to keep the 
entire e expansion in a single place. This procedure is even more motivated when we 
absorb these extra terms into L\ and L 2 , for these terms are necessary in order to 
satisfy L-\U\ = 0. If instead we kept separate the e expansion of u \, the result would be 
an awkward expression like L 0 -u 1;1 + LiUi- 0 = 0. To introduce our convenient shortcut, 
first recall that since N\ is never needed, we only need to perform this within L. Then, 
for the places where U(ij ± |) occurs within L, we substitute, keeping only to the 
order required, 


y 2 edl^j 

and then later on we substitute the specific form of U\, we substitute A(Y) rather 
than A(Y + \A^ 2 y). The second place we introduce the expansion is that since r 
depends on dy, we have 

f) 1 f) 2 

r{d ¥ ) —► r(dy + e 1/2 d Y ) = r(dy) + e 1/2 d Y -^^r(dy) + - e^-^^r(%). (C.17) 

Then, letting dy —>■ iq (which we can do because we will only need to perform this 
expansion on a term e iq, -' y and not other harmonics), we see that 

/ d ^ d 2 \ 

r(q) -)■ ^1 - ie ll2 d Y — - -ed^-^jr{q). (C.18) 

To implement this, one can set, in L, 

W(y) ^ (l - ie l,2 d Y “d q ” - ^ed$ u d 2 ”^W(y), (C.19) 

where the “ d q v means that the derivative acts only on r(q), not on the e iqy part of 
U\ . We need only make this replacement in W, not U, because the U component of 
the right null eigenvector does not contain any derivatives dy. One can verify that 
this shortcut gives the same results as if one proceeded more straightforwardly. 

The problem is most conveniently expressed in terms of the Fourier transform of 
the difference variables x, y. We use the convention 

f{k x ,ky) = I dxdye- ik * x e- lk yyf(x,y). (C.20) 


U(y±\y), (C.16) 


U(y± \y) ->• (l ± ^ ye 1/2 d Y + ^ 
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After Fourier transform, the required linear operators are given by 


[L 0 u\w — 
[L 0 u]u = 


{ 2/3k x kydy 
V 9o(k x ,ky,y) 



W(k x , k y | y) 


dk y e 


ik : 


' yV h 0 (k x , k y , k 1 


2vr 


d„ 


f dKdk y —^l 

(2vr yj go\k x , k y , y) 


--W(k x ,k y | y) - yU{y), 


(C.21a) 

(C.21b) 


[L 2 u\ w = -d T W(k x ,k y | y) - dk' e lk 'y y h 0 (k x , k y , k 


,AK) 


2tc 


‘Ifidy ( Soifoxi ky-> y) 


ik x k y d‘yQ\ (k x , k y , y) 




2.9o(fc 3 


[L, 2 u]u — —drU(y) + d\ 


go(k x , ky, y) 

&-j } ‘‘d^W(k„k y I y)+c^y“Bf’W(K,k, I y), (C.22a) 

sofoxj ky, y) 


(2t) 


dk x dk y 


ik x kydygi(k x , k y , y) 


“ d q ” + 


v x^yvy u q2v 


go(k x ,k y ,y) q 2g 0 (k x , k y , y) q \ 
And the nonlinear operator is given by 


W(k x ,k y \y). (C.22b) 


[iV 0 un (u,z)] w = J dky e lk y y si(k x , k y , k' y , y)z w (k x , k y - \k' y \ y) 

Vu(ky. 


- Si(k x , ky , -k' y , —y)z w (k x , k y + \k' y \ y) 


2tt 


[Nr(v,z)]u = 0, 


(C.23a) 

(C.23b) 


where again the symmetrized version is N 0 (v,z) = ^[N^ n (v, z) + Nq u (z,v)}. 

Here, U is the Fourier transform of U. Our expressions for U(y) will always 
consist of periodic exponentials, and so U contains delta functions and the convolution 
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integral can be immediately performed. We also have defined 


k 2 = k 2 x + k 2 , 

9o{k x , k y , y) = k A - ^kid 2 + i k 2 yd 2 + d 
Qiikxi k y , y) = -k 2 x + k 2 y + 


9 - 2 (k x ,k y ,y) = ~k 2 x + \k 2 y + \d 2 


ho(k x , ky, ky) 


k' 2 

Ky 


4 V 

2 n y + 8 ^’ 

k\) ( k X , ky 2 ky ) 


(C.24) 

(C.25) 

(C.26) 

(C.27) 


k 2 x + (k y - \k' y ) 2 

F 0 (k x ,k y + \k ' y ) }•, 


U2 

K y 


k 2 x + (k y + \k 


1 U\2 


(C.28) 


So(^a;5 ky, y ) k x ky8y 


dy9\(k x , ky, y) 2 - g 0 (k x , k y , y)g 2 (k x , k y , y) 9l (k x , k y , y) 


9o(k x ,k y ,y ) 3 


g 0 (k x ,ky,y) 2 

(C.29) 

. ,f. h, *?[-*£ - (*» - h;) 2 +i 9 ?] --1 4 ;) 

(C.30) 

We also define go(k x , k y , ky) as go(k x , k y , y) with <3y —» ik y , and similarly for gi, g 2 , sq, 
and S i. 

We dehne an inner product 

(v,z)= dyvu(y)zu(y) + / dy dk x dk y v* w (k x , k y , y)z w (k x , k y , y). (C.31) 


At 0(U/ 2 ), we hnd that U\ is given by (C.6) where the right null eigenvector is 
given by r = e iqcV {r w (k x , k y , q), Uq), where 


rw(k x , k y , q) 


h(k x ,k y ,q) 


hp(k x , k y , q)Up 
g 3 (k x ,k y ,q) 

+ 2 i/3k x k y q 
9o(k x , k y , q) ’ 


(C.32) 

(C.33) 


and Uq is a constant with dimension of velocity, whose purpose is to help keep track 
of dimensional consistency. For any computation it can be set to unity. The complex 
conjugate of the left null eigenvector is found to be l* = e~ tqcV {l) v (k x , k y ), 1}, where 


7 * (k k) = _ ^kxky _ 

w *’ v (27T) 2 g 0 (k x ,k y ,q)g 3 (k x ,ky,q) 


(C.34) 
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The q dependence of rw and lw is now suppressed except for where it matters in 
(C.45); they should be evaluated at q — q c . 

At 0(e), we need to solve the particular solution of (C.9). Take an ansatz 


u 2 pW — a w(k x , k y ) A 2 e i2qcy + aw(k x , k y )*A* 2 e bw(k x , k y )AA*, (C.35a) 

u 2pU = a u A 2 e l2q, ' v + a^V* 2 ** + b v AA*. (C.35b) 

After some algebra we find 

/ 7 r \ _ h 0 (k x , k y , 2q c ) 

aw(k ” K) ~ ~HK,k v , 

+ U 0 g 3 (k x , k y , 2 g ) _1 [§i(k xi ky, Qci q c )r w (k XI ky 2 ^ C ) 

- Si(k x , k y , -q c , —q c )r w (k x , k y + \q c )\, (C.36) 

Na 
D a ’ 

2 iq c U 0 


au = 


(C.37) 


N = 
a (2tt) 2 


dk x dk y k y , 2q c )gs(k x , k y , 2 (^ C )J 


-i 


x [Si(A:*, k y , q c , q c )r w (k x , k y - \q c ) - s x (k x , k y , -q c , -q c )r w (k x , k y + \q c )\, 

(C.38) 


n _ , 2 0 , 


k x k y h 0 (k x , k y , 2q c ) 


c dk dk - 

x y g 0 (k x ,k y ,2q c )g 3 (k x ,k y ,2q c y 


(C.39) 


and 


On 


b w {k x , k y ) = ^ [hiK, k y , q c , -q c )r w (k x , k y - \q c )* 

- h{k x , k y , —q c , q c )r w (k x , k y + \q c )* - h{k x , k y , q c , -q c )r w {k x , k y - \q c ) 


+ si (k x , k y , —q c , q c )r w {k x , k y + \q c )], 
bu = 0 . 

At 0(e 3//2 ), the solvability condition (C. 10 ) becomes 

c 0 d T A(y, t) = ci A + c 2 dyA - c 3 |A| 2 A, 


(C.40) 

(C.41) 


(C.42) 
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where 


Co Uq T I dk x dk y l w {k x , k y )rw{kx, ^/) | g=?c > 


Ci Uq I dk x dky l w (k x , k y )ho(k x , fcy) | ?=g , 


c 2 = 


C 3 = 


<9g 2 


I dk x dky l w (k x , k y , q)Uoh 0 (k x , k y , q), 


(C.43) 

(C.44) 

(C.45) 


q=q c 


dk x dky l\y(k x , ky ) s Uq |^si (k x , ky , q, 0 ) 6 jy {kxi ky 2 *?) 


-si(/c x ,fcy,-g, 0)6 h/(/c x ,^ + !<?)] +f/ 0 [si(/c a; ,fcy,-g,2g)a H /(/c a; ,A:y + |g) 
- Si(A; x , k y ,q, -2 q)a w (k x , k y - \q)\ + [s 1 (fc a: , k y ,2q, -q)r w (k x , k y - q)* 


- h{k x , ky, -2q, q)r w {k x , fcy + g)*] | 


(C.46) 


<?=9c 


After returning to the unsealed variables by letting T —y et, Y —y and 

A —y A/e, we recover (3.60). The coefficients c t involve integrals over the forcing 
spectrum which is here presented in a form where the the wavenumbers are shifted, 
i.e., contain terms like k y — q/2. It is also possible to shift the integration variable so 
all integrals contain just the unshifted forcing F 0 {k x , k y ), after which the q derivatives 
in C 2 can be explicitly computed (Bakas and Ioannou). 

It is possible to obtain Co, Ci, and C 2 , which govern the linear behavior, via the 
alternate and much simpler route of using the analytic dispersion relation (3.25). 
The dispersion relation can be put into the form D( A, e, q) = 0. The conditions of the 
instability threshold require D( 0, 0, q c ) = 0 and dD/dq( 0, 0, q c ) = 0. Thus, expanding 
D to lowest order about (0, 0, q c ), we find 


dD 

~d\ 


(0, 0, q c ) A 



(0,0, q c ) + 


1 d 2 D 

2 ~d(f 


(0, 0, q c )(q 



(C.47) 


Then up to a constant of proportionality, we see that c 0 = —dD/d\(0,0,q c ), C\ = 
dD/de(0,0,q c ), and c 2 = — /d 2 D/dq 2 (0, 0, q c ). This was used to put c 2 above into 
a succinct form. But this approach does not give c 3 ; for that one needs the full 
bifurcation calculation which includes nonlinear terms. 

To verify these analytic expressions, we take an example forcing F 0 (/e x , k y ) = 
Ak exp [—(k — kf) 2 /al], with k 2 = k 2 + k 2 , A = Ay/ne/ak, kf = 1, and cr k = 0.5. For 
the other parameters we use n = 0.1, /3 — 1. Then the critical value of the control 
parameter is calculated to be e c = 0.1297 with a critical wavenumber q c = 0.676. 
We compute cq = 1.10, ci = 0.10, C 2 = 0.0015, and c 3 = 2.28. Comparisons between 
analytic and numerically computed results are shown in Figure 3.8 and are in excellent 
agreement. 
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Appendix D 

Projection for Ideal State 
Equilibrium 


In this Appendix we perform the projection of (4.1) onto the Galerkin basis functions. 
We find explicit formulas for the nonlinear algebraic equation, in a suitable form 
for numerical implementation. The shorthand notations in (4.21)-(4.29) are used 
throughout. 


D.l Eddy Equation 

Projection of the eddy equation entails operating on (4.1a) with 


27r27r2vr\" 1 r /a r /b r /q _ . 

-,-/ dx dy dy <j> rat , 

a b q J J—Tr/a J~ir/b J—ir/q 


where 


■'mnp 


— gimax e inby £ ipqy 


(D.l) 


(D.2) 


First Term 

First term of (4.1a): 

~[U+ ~ U-}d x W = -[U(y + \y) - U(y-\y)]d x W. (D.3) 


Note that 


p 

U(y + \y) - U(y - \y) = ^ U p ,e ip ' q y {e ip ' qy/2 - e ~ ip ' qvl2 ) (D.4) 

P '=-p 


and 

d x W = Y, W mn P 4e“e^e«. (D.5) 

mnp 
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Then we have 


■[U+ - U-]d x W = -J2 U p 'W mnp ik x e ip ' qy (e ip ' qy/2 - e - ip ' qy / 2 )e imax e inby e ipqy . (D.6) 


p'mnp 

Now project onto 4> rstl yielding 


U T yW mnp lk x ^ [ n/a dxe^-r^-L- f' b dye Mnb- S b)^ P ’ qy /2 _ g— ip'i} 2 // 2 ) 


p'mnp 


: 2n/a J_ n/a 


2 tt/6 J_ n/b 

^ j-ir/q 


X 


'in Ji.p'+p-tfqy 


2 tt /q J_, 


/<? 


dy e 


For the y integral, use 


pn/b 


/ dy e iay = sine | 

A 

/ — 7 r/b 


Performing the integrals results in 


with 


rstp'mnp 

Now split into real and imaginary parts: 


t(X) TT T/T/ 
1 rstp'mnp u P l vv ™np, 


1rstp'mnp ~ ik x ((J+ (J-)8 mir 8 p ' + p-t, 0 - 


where 


r(l) TT TT/ _ 7(1) TT 77 1-77(1) TT Tp 

1 rstp'mnp u P ,yy mrip ~ J rstp'mnp u P 1 r ™rip T bl ^ r stp'mnp u P 1 ^mnp, 


^rstp'mnp ^rstp'mnp k x [(J-\- tfi 


Second Term 

Second term of (4.1a): 


l4-l/':)(v 2 +is§ ]£>,!■. 


where 


And 


U" ± =Ul- a ZF L~ 2 U ± = -J2 U A 


’ 2 pip'qy P ±w'qy /2 
y,u e e 


^ + -d 2 )d x q = -J2 C ™n P ik x [k + -kl w 


mnp 


t-2 1 . 


£ imax inby £ ipqy 


(D.7) 


(D.8) 

(D.9) 

(D.10) 

(DU) 

(D.12) 


(D-13) 

(D.14) 

(D.15) 
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Then 


(V+ -f/")(v 2 + \e%\djs = Y u p/ c mnp ik x klJk 2 + 

\ / p'mnp ' ' 

x ^ e ip'qy/2 _ g—®P , 92//2^ e imax £ inby £ ipqy 

Project onto (p rst and obtain 


d 2 ) TT (~f 

I rstp'mnp u p'^mnpi 


with 


Irstp'mnp — ikxky u + ^/Cyjy^((T + CT_)<5 mjr 5p'+p_t,l 


Now split into real and imaginary parts: 


r( 2 ) TTr 1 — A J ‘ > TT TJ i zyOO tt rr 

1 rstp'mnp u p'^ y mnp <J r stp'mnp u p’ 11 mnp ~r l ’ 1 '-rstp'mnp' J p'^mnp-, 


A 2 ) 


'(2) 


(D.16) 

(D.17) 

(D.18) 

(D.19) 


where 


t( 2 ) _ _ 7^(2) _ 

u rstp'mnp v rstp'mnp 


= -kxk 2 y tU (V + (<T+ - (T_)5 mir 5p/ +p _ t)0 . (D.20) 


Third Term 

Third term of (4.1a): 

2/3dyd y d x ^. 

The projection process should be clear. Here we obtain 

r(3) 

1 rstmnp^'map i 


with 


1rstmnp T‘2(3k x kyky j w8 rn r8 n: s8p : t 


In real and imaginary parts: 


r( 3 ) _ tW tt i • t>WU ^ 

± rstmnp K -mnp °rstmnp 11 mnp t ul '-rstmnp yJ mnp 


7(3) 


'(3) 


with 


^rstmnp Krstmnp ‘^ , P^x^y^y,W^m,r^n,s^p,t 


43) 


Fourth Term 

Fourth term of (4.1a): 

- (u+ + U"^jdyd y d x ^>. 
After projection onto (f> rs t , we obtain 

T ( 4 ) tt ri 

1 rstp'mnp u p' ^mnpi 


(D.21) 

(D.22) 

(D.23) 

(D.24) 

(D.25) 

(D.26) 

(D.27) 
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with 


r(h 


J rstp'mnp 

Now split into real and imaginary parts: 


^ky U^x^y^y,W^+ Y ^—)^m,r^p'+p—t, 0- 


r(4) TT C 1 — T^ TIN zvd 4 ) TT C 1 

1 rstp'mnp' J P , ^'mnp J r stp'mnp u p ' 11 mnp T rstp'mnp u p'^mnpi 


where 


7 ( 4 ) 

^ rstp'mnp 


—K^ 4) 

rstp'mnp 


i~2 


^y,U^ x ^y^yiW(^+ “I - &—^)$m,r3p'-\-p—t, 0* 


(D.28) 


(D.29) 


(D.30) 


Fifth Term 

Fifth term of (4.1a): 

F(x,y). (D.31) 

This is the forcing term. After projection onto 4> rs t, we obtain 

42 = 42 = <5t,o (y y) F nb(k x , ky). (D.32) 

where F nh is “narrowband forcing” as described by Srinivasan and Young (2012) in the 
continuous Fourier transform. The prefactor of {2tt / o)~ 1 {2tt/ b)~ x essentially comes 
from the conversion factor from a Fourier transform amplitude to a Fourier series 
amplitude. Specifying F n y(k Xl k y ) defines F(x,y). We use 


Fnb(k X i ky ) 


2'Kekf/8k kf — 5k < k < kf + 8k 
0 otherwise 


(D.33) 


where £ is an equivalent energy input in the case of L d 2 = 0. Note is pure real. 


Sixth Term 

Sixth and final term of (4.1a): 


— (2/r + 2uDy)W. 


After projection onto (j) rs t we obtain 


7 ( 6 ) 

1 rstmnp 


IT 


mnp? 


with 

42 mnp — — [2/r + *4 h ^ 

In real and imaginary parts: 


7(6) 

1 rstmnp 


IT, 


mnp 


T ( 6 ) pr 1 j zvw F 

° rstmnp -‘-'mnp ~ rstmnp - 1 mnpi 


46) 


(D.34) 


(D.35) 

(D.36) 

(D.37) 
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with 


Jrttmnp ~ ^rltmnp ~ ~ ^A 4 + + h 2 -)]^m, r ^n,s^p,t- (D.38) 

D.2 Zonal Flow Equation 

For the zonal flow equation (4.1b), we can project onto e ltq or equivalently, just 
equate the e ipgp coefficients. Upon substitution of the Galerkin series, we have 

-U p (ji + uk^k'y jj/k^ u + ^ ik x k y kyC mnp = 0, (D.39) 

mn 

where here we use kyjj = pq. Noting the real part of C cancels out after summation, 
this becomes 

-Up(/i + vkf)k 2 y,u/kyjj ~ ^ k xk y kyH mnp = 0. (D.40) 

mn 

By symmetry, one could sum only over positive m, n (and put in a factor of 4), though 
we do not need to do this. The above equation can be written 

Itp'Up' + ItmnpHmnp = 0, (D.41) 

where 

I?J = -(/* + vl$u) {klu/klu)^P', (D.42) 

1 tmnp k X kyky\y5t p . (D.43) 
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Appendix E 

Projection for Ideal State Stability 


In this Appendix we perform the projection of the linearized system (4.53) onto the 
basis functions. We find explicit formulas for the matrix equation, in a suitable form 
for numerical implementation. The shorthand notations in (4.66)-(4.79) are used 
throughout. 

E.l Eddy Equation 


Let 


- / mnp 


— £ imax £ inby £ i(Q+pq)y 


(E.l) 


We will project (4.53a) onto 4> rst in the same way as for the ideal state calculation, 
by applying 


/ 2tt 2tt 27t\ -1 
\ a b q J 




dy 0*rst 


(E.2) 


We suppress the e at dependence of the perturbations 5W and 5U from now on. Pro¬ 
jecting the LHS of (4.53a) is trivial; one merely obtains crSW rst . Now we project the 
RHS. The matrix coefficients are closely related to those in the ideal state calculation. 


First and second term 

The first term on the RHS is 


-{SU+ - SU_)d x W. 


(E.3) 


We have 


so that 


5U ± = 5U(y±\y)=Y J dU P >e i 


(Q+p'q)y e ±i(Q+p'q)y/2 


5U+ — 5U- = ^ SU f 


/e 


i{Q+p’q)y (J{Q+p'q)y/2 _ -i{Q+p'q)y/2 


(E.4) 

(E.5) 
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Also, 


d x W = Y, Wmn P ik x e imax e inby e ipqy . (E.6) 

mnp 

Thus, 


-(SU+ - 5U_)d x W = - SU p W mnp ik x e^ Q+p ' q)y e imax e inby e ipqy 

p'mnp 

x ^ e 7 < 2+r / '?)y/ 2 — e -i(Q+p'q)y/2' s j _ (E.7) 

Now project onto (p rst . Obtain 

I% m n/U^W mnp , (E.8) 

where 

Irstp’mnp ~ ‘~‘ik x ((J+ t 8U ~ cr -,Su)^m,r^p'+p-t,0- (^-9) 

The second term on the RHS of (4.53a) is 

-(17+ - UJ)dJW- (E.10) 

With the aid of our convenient notation, we can obtain the result after projection 
from the first term by making the replacements 5U —> U and W —> SW (including in 
the subscripts of the coefficients). We obtain 

m,nn p U^W mv , (E. 11) 

where 

~/ 9 \ 

Kstp'mnp = —ik x (cr +t u — cr_ ! [/)5 mir (5p'+p_t i 0. (E.12) 


Third and Fourth Term 

Third term: 

(sul - tu "_) (V + d x y. (E.13) 

We obtain after projection 

/U^ mnr , (E. 14) 

where 

^rstp'mnp ~ ^xky^gu ^k + ~^k p W ^j (0+,<S£/ — G- i su)^m,r^p'+p—t,0- (E.15) 

The fourth term is obtained after the appropriate replacements: 

^rstp'mnpUp'^mnp, (E.16) 
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where 


Kstp'mnp i'kxky uyk + ky {5W ) (&+,U cr —) < ^r7T,,r' < 5^p / +p_i,0- 


Fifth and Sixth Term 

Fifth term: 

We obtain after projection 
where 


8U++8U"_ )dyd y d x q. 


7t5) rrr vjy 

1 rstp'mnp ^ P' mnpi 


Irstp'mnp ^y,Sukxkyk yj w{^+,SU + cr -,5U') ( 5m,r ( 5p'+p-t,0- 

The sixth term is obtained after the appropriate replacements: 

m 


.-r-2 


T(°) TT rrii 

1 rstp / mnp u P lu ^mnp, 


where 


I rstp'mnp ~ * ky U k x k y ky } Sw{ (J +,U + & - ,u)^m,r^p' +p-t,0- 


—2 


Seventh and Eighth Term 

Seventh term: 


(E.17) 


(E.18) 

(E.19) 

(E.20) 

(E.21) 

(E.22) 



2/3dyd y d x S'i>. 

(E.23) 

After projection, we obtain 

?( 7 ) Avtf 

1 rst.mnp u 1 mnpi 

(E.24) 

where 

1 rstmnp ^P^x^y^y,5W^m,r^n,s^p,t' 

(E.25) 

Eighth term: 

— (2 Jl + 2 uDh)SW. 

(E.26) 

After projection we obtain 

7 8 ) AW 

1 rstmnp u vv mnpi 

(E.27) 

where 

7(8) 

1 rstmnp 

- [2 11 + l/(h^g W + 

(E.28) 

E.2 Zonal Flow Equation 


Since the zonal flow equation (4.53b) is linear, projection is equivalent to matching 
the coefficients of the exponentials. It is simple to fold that at each p, 

h 2 

XTT ■ fu y, 8U 

aoU p = i ^— 

ky,5U 

^ ^ kxkyk'y,mnp (/^ “f - ^ / ^ , y i 8u)^^P’> 

mn 

(E.29) 
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_ 2 

where here we use ky t $u = ky,sw = Q + pq , k x = ma, k y = nb, and ky 5u = k 
&zfL^ 2 . 
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Appendix F 


Stability of Kraichnan-Spiegel 
Closure 

F.l The Kraichnan-Spiegel Model 

In the pedagogical closure described in Chapter 5, we proved that any nonzero steady- 
state solution is linearly stable. Here, we prove the same result for the Kraichnan- 
Spiegel (KS) closure (Kraichnan and Spiegel 1962). The KS closure is a model for 
energy transfer in 3D, isotropic turbulence. The significance of this result stems from 
the fact that stability of the equilibria of closures is an important topic, and this is 
one of the first definitive results. 

A limit of the KS closure assuming local transfer gives the Leith diffusion model 
(Leith 1967). Hence, stability of the Leith model follows from stability of the KS 
model. 

In general, an energy balance equation can be written 

= 2 lt E(k) + T(k), (F.l) 

where E(k ) is the energy spectrum and T(k) is the nonlinear transfer term. The 7 *, 
term includes all linear terms, generalizing the original KS model by allowing not 
only for viscous damping but linear drive as well. In Navier-Stokes, the quadratic 
nonlinearity means that in fc-space the fundamental interactions are among three 
Fourier modes, or triads. The KS approximation involves treating the fundamental 
nonlinear transfer as occurring only between two modes. The KS closure takes a 
specific form for T(k): 

dE g^ = 27 kE(k) + dp [S e (p I k) - S e (k I p )], (F.2) 

where the “emission” term S e (k \ p), corresponding to the energy emitted by mode k 
and absorbed by mode p, is given by 

S e {k | p) = V [k- 2 E(k)]- 3/2 (kpY/ 4 g(p/k), 
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(F.3) 




and S e (p | k) is an absorption term at mode k corresponding to emission from mode 
p. Here, rj is a dimensionless numerical constant, g(p/k ) = g(k/p), and g(x) decays 
quickly for x 1 (g enforces locality in wavenumber space). We will not be concerned 
with the functional form of g(k/p), only its symmetry, so we write g{k/p) —> g kp and 
note that it is symmetric in its indices. The balance equation written explicitly is 

= 2 7 k E(k)+ V J dpg kp (kp) 7 / A [p- 3 E(p) 3 / 2 - k~ 3 E(k) 3 / 2 ]. (F.4) 

Now, assume that a steady state solution E{k) exists which is nowhere zero. Then 
one can write a “steady-state condition” which will be later used to eliminate 7 k : 

2lk = W)J dv9kp {kp)7/4 \ k ~' E ^' 2 - p~ 3 e (p) 3/2 } ■ ( F - 5 ) 

F.2 Linear Stability 

We show that in the Kraichnan-Spiegel closure, any nonzero equilibrium is linearly 
stable. We do this by providing a positive definite functional, quadratic in the per¬ 
turbation, which decays in time. 

Linearization 

Linearize about an equilibrium E(k), assuming one exists, by letting E(k , t) = E(k ) + 
5E(k,t). Then (F.4) becomes 

^ = 27 kSE(k,t) + y J dpg kp ( kp) 7/ 4 

x [p- 3 E(p)^ 2 SE(p, t ) - k~ 3 E(k) l ' 2 8E(k 1 1)]. (F. 6 ) 

It will be convenient to write this in terms of w k = 5E(k,t)/E(k): 

= 27 k w k + J dp g kp ( kp) 7/i [p~ 3 E(p) 3/2 w p - k~ 3 E{k) 3/2 w k ]. (F.7) 

Now, substitute the steady-state condition (F.5) to obtain 

aj w = W)S dpgkr (kp}7/ ‘ [*" 3B w3/2 ^ p' 3e ^ 3 ' 2 ] w ‘ 

+ 2E(k) / d ' t,9t r( k P) lli [P~ 3E ('l > f l2w p~ k~ 3 E(kf l2 w k }. (F. 8 ) 

The notation can be simplified by defining 

Q k = k~ 3 E(k) 3 / 2 , 

Ck p = g kp ( kp ) 7/4 , 
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(F.9) 

(F.10) 










where C kp is symmetric in its indices, to yield 


dw k 

dt 


V 

2 E(k) 


J dp C k p (^Q k Wk T ^QpWk 3Q p Wp'j. 


(F.ll) 


Quadratic Functional 

Consider the quadratic functional 

1 f°° 

W{t) = ^J dk Q k E(k) w\. (F.12) 

XV (t ) is positive definite with respect to w k - The evolution of XV {t ) is given by 


dXV(t) 

dt 


J dkQ k E(k) w k^^~ 

2 dP dp E k p ( Q k ^ k “I - 2QkQ P w k ^QkQp^k^p) • 


(F.13) 


We will now show that dXV/dt < 0, meaning that perturbations decay and the 
equilibrium is linearly stable. Note that for terms inside the square brackets in (F.13), 
we are free to swap the indices k v-)- p, since C kp is symmetric in k,p. Then, through 
a series of manipulations using this fact (we use the equals sign as if the following 
took place under the integral), 

QkW k T 2QkQpW k 3Q k QpW k Wp 

QkW k T QkQpW k “I - QkQptVp SQkQpWkWp 

QkQpW k ^QkQpW k Wp T QkQpWp T Qkd^k Q kQ pW k Wp 

QkQp(u>k 2tCfctUp d" Wp) T Qk^k QkQpWkWp 

QkQp(w k w.p) T Qk^k QkQpWkWp. (F.14) 

In the second line, we have let 2 Q k Q P w k QkQ P w\ + QkQpWp. Now, observe 

Qk^k QkQpWkWp 2 QkW k ^QkQpWkWp) 

2 {.Qk^k 2 QkQpW k Wp d" QpWp^ 

= | (QkW k - QpWpf. (F.15) 


Putting this all together, we obtain 

QkQpiWk Wp) T — i^QkWk QpWp) 

We have found that dXV/dt < 0, and vanishes only when w k = 
required any further conditions on g(x). 


dw n f 
— = / iki P C kr 


< 0. (F.16) 

0. We have not 
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With external forcing 

It is not difficult to add random (isotropic) forcing to the model, which becomes 
nonrandom, positive F(k) at the energy balance equation. The conclusion remains 
unchanged, as we now show. The energy balance equation can be written 

=F(k) + 2-t k E(k)+v [■■■ , (F.17) 


where the • • ■ indicate terms that were previously present without forcing. The steady 
state condition, assuming E(k) is nonzero everywhere, is now 


2 q k = 


F(k ) 


+ 


E(k) ' E(k) 

In the linearization equation, the F(k) vanishes, giving 


dw k 

dt 


= 2~f k w k + 


3 rj 


2 E(k) 


Substituting the steady state condition gives 

dw k _ F(k) _ rj r 

dt ~ E(k) Wk 2 E(k) J 


(F.18) 


(F.19) 


(F.20) 


Taking the same quadratic functional W(t), we obtain the evolution equation 

dw r 

_ = _y dkF(k)Q k wl + --- . (F. 21 ) 


The new term involves the forcing F(k). In dW/dt , the forcing contributes a negative 
definite term, so the total dW/dt is still negative definite. 
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